Fair sampling of ground-state configurations of binary optimization problems

TL;DR

This paper introduces a novel sampling method combining parallel tempering Monte Carlo and isoenergetic cluster moves to fairly sample ground states in binary optimization problems, addressing biases in existing heuristics.

Contribution

The paper presents a new approach for unbiased sampling of degenerate ground states in binary optimization, applicable to quantum annealing and classical heuristics.

Findings

Method effectively samples degenerate states with nearly equal probability.

Applicable to Ising spin glasses and quantum annealers like D-Wave.

Includes a heuristic for estimating the number of solutions.

Abstract

Although many efficient heuristics have been developed to solve binary optimization problems, these typically produce correlated solutions for degenerate problems. Most notably, transverse-field quantum annealing - the heuristic employed in current commercially-available quantum annealing machines - has been shown to often be exponentially biased when sampling the solution space. Here we present an approach to sample ground-state (or low-energy) configurations for binary optimization problems. The method samples degenerate states with almost equal probability and is based on a combination of parallel tempering Monte Carlo with isoenergetic cluster moves. We illustrate the approach using two-dimensional Ising spin glasses, as well as spin glasses on the D-Wave Systems Inc. quantum annealer chimera topology. In addition, a simple heuristic to approximate the number of solutions of a degenerate problem is introduced.

Tables

Table 1. Table 1: Parameters of the simulation: For each instance class and system size N 𝑁 N , we compute N sa subscript 𝑁 sa N_{\rm sa} instances. N sw = 2 b subscript 𝑁 sw superscript 2 𝑏 N_{\rm sw}=2^{b} is the total number of Monte Carlo sweeps for each of the 4 ​ N T 4 subscript 𝑁 𝑇 4N_{T} replicas for a single instance, T min subscript 𝑇 min T_{\rm min} [ T max subscript 𝑇 max T_{\rm max} ] is the lowest [highest] temperature simulated, and N T subscript 𝑁 𝑇 N_{T} and N hc subscript 𝑁 hc N_{\rm hc} are the number of temperatures used in the parallel tempering method and in the isoenergetic cluster algorithm, respectively.

Topology	Couplers	$N$	$N_{\mathrm{sa}}$	$b$	$T_{\min }$	$T_{\max }$	$N_T$	$N_{hc}$
2D	$\{\pm 1,\pm 2,\pm 4\}$	$144$	$360$	$24$	$0.05$	$3.05$	$35$	$35$
2D	$\{\pm 1,\pm 2,\pm 4\}$	$256$	$360$	$24$	$0.05$	$3.05$	$35$	$35$
2D	$\{\pm 1,\pm 2,\pm 4\}$	$576$	$322$	$24$	$0.05$	$3.05$	$35$	$35$
2D	$\{\pm 1,\pm 2,\pm 4\}$	$784$	$232$	$24$	$0.05$	$3.05$	$35$	$35$
2D	$\{\pm 1,\pm 2,\pm 4\}$	$1024$	$370$	$24$	$0.05$	$3.05$	$35$	$35$
Chimera	$\{\pm 1,\pm 2,\pm 4\}$	$128$	$360$	$24$	$0.05$	$3.05$	$35$	$20$
Chimera	$\{\pm 1,\pm 2,\pm 4\}$	$288$	$360$	$24$	$0.05$	$3.05$	$35$	$20$
Chimera	$\{\pm 1,\pm 2,\pm 4\}$	$512$	$360$	$24$	$0.05$	$3.05$	$35$	$20$
Chimera	$\{\pm 1,\pm 2,\pm 4\}$	$800$	$360$	$24$	$0.05$	$3.05$	$35$	$20$
Chimera	$\{\pm 5,\pm 6,\pm 7\}$	$800$	$976$	$24$	$0.10$	$1.55$	$30$	$23$
Chimera	$\{\pm 1,\pm 2,\pm 4\}$	$1152$	$223$	$24$	$0.05$	$3.05$	$35$	$20$

Table 2. Table 2: For each instance with system size N = 1024 𝑁 1024 N=1024 , we run N sw = 2 23 subscript 𝑁 sw superscript 2 23 N_{\rm sw}=2^{23} Monte Carlo sweeps for each of the 4 ​ N T = 4 ​ N h ​ c = 120 4 subscript 𝑁 𝑇 4 subscript 𝑁 ℎ 𝑐 120 4N_{T}=4N_{hc}=120 replicas with lowest temperature T min = 0.17 subscript 𝑇 min 0.17 T_{\rm min}=0.17 and highest temperature T max = 1.3 subscript 𝑇 max 1.3 T_{\rm max}=1.3 . The fixed number of Monte Carlo sweeps and free variables for each iteration are 1 / 24 ​ N sw 1 24 subscript 𝑁 sw 1/24N_{\rm sw} and 1 / 24 ​ N 1 24 𝑁 1/24N , respectively. The median estimate of the degeneracy G RG superscript 𝐺 RG G^{\rm RG} is averaged over 10 10 10 independent runs and error bars are computed using the jackknife method.

Instance	$G^{\mathrm{exact}}$	$G^{\mathrm{RG}}$	error	% error
$0$	$2.0094\times {10}^{29}$	$1.9415\times {10}^{29}$	$\pm 6.04\times {10}^{27}$	$3.00\%$
$1$	$9.7771\times {10}^{34}$	$1.0081\times {10}^{35}$	$\pm 3.58\times {10}^{33}$	$3.66\%$
$2$	$3.3778\times {10}^{27}$	$3.3188\times {10}^{27}$	$\pm 1.18\times {10}^{26}$	$3.50\%$
$3$	$1.2826\times {10}^{32}$	$1.3041\times {10}^{32}$	$\pm 2.57\times {10}^{30}$	$2.00\%$
$4$	$1.8613\times {10}^{39}$	$1.9317\times {10}^{39}$	$\pm 7.59\times {10}^{37}$	$4.08\%$
$5$	$1.4104\times {10}^{40}$	$1.4515\times {10}^{40}$	$\pm 7.28\times {10}^{38}$	$5.16\%$
$6$	$9.6510\times {10}^{29}$	$9.6105\times {10}^{29}$	$\pm 1.83\times {10}^{28}$	$1.90\%$
$7$	$2.3699\times {10}^{38}$	$2.3543\times {10}^{38}$	$\pm 1.66\times {10}^{37}$	$7.04\%$
$8$	$1.4168\times {10}^{31}$	$1.3527\times {10}^{31}$	$\pm 1.16\times {10}^{30}$	$8.58\%$
$9$	$1.1265\times {10}^{34}$	$1.0789\times {10}^{34}$	$\pm 5.06\times {10}^{32}$	$4.69\%$

Equations

$$\mathcal{H}=-\sum _{i<j}^N{J}_{ij}{s}_i{s}_j,$$

$$Q_{\mathrm{th}}=\sigma /e=\sqrt{(1-p)/np}=\sqrt{(G-1)/n}.$$

$$R_{i-1}=G_{i-1}/G_i$$

$$G^{\mathrm{RG}}=G_{\mathrm{final}}\prod _i{R}_{i-1}.$$

Full text

I Introduction

Quantum annealing Finnila et al. (1994); Kadowaki and Nishimori (1998); Brooke et al. (1999); Farhi et al. (2001); Santoro et al. (2002); Das and Chakrabarti (2005); Santoro and Tosatti (2006); Das and Chakrabarti (2008); Morita and Nishimori (2008) and, in particular, quantum annealing machines have ignited an ever-increasing interest in algorithms used in statistical physics to solve hard combinatorial industrial optimization problems, as well as related applications. While there has been an extensive body of work attempting to discern if the D-Wave Systems special-purpose quantum annealing machine can outperform algorithms on conventional CMOS hardware Dickson et al. (2013); Pudenz et al. (2014); Smith and Smolin (2013); Boixo et al. (2013); Albash et al. (2015a); Rønnow et al. (2014); Katzgraber et al. (2014); Lanting et al. (2014); Santra et al. (2014); Shin et al. (2014); Boixo et al. (2014); Albash et al. (2015b, a); Katzgraber et al. (2015); Martin-Mayor and Hen (2015); Pudenz et al. (2015); Hen et al. (2015); Venturelli et al. (2015); Vinci et al. (2015); Zhu et al. (2016); Mandrà et al. (2016); Mandrà and Katzgraber (2017); Mandrá and Katzgraber (2018), there have been only a few studies Boixo et al. (2013); Albash et al. (2014); King et al. (2016); Mandrà et al. (2017); Könz et al. (2018) attempting to characterize the sampling ability of quantum annealing. Initial studies Matsuda et al. (2009); Mandrà et al. (2017) suggested that transverse-field quantum annealing with stoquastic drivers result in biased solution distributions for degenerate problems. However, more recently, it was shown Könz et al. (2018) that even with high-order drivers the sampling bias can be removed only in special cases.

Many industrial applications rely more on a broad solution pool then on the minimum of the cost function, with some prominent examples being propositional model counting and related problems Jerrum et al. (1986); Gomes et al. (2008); Gopalan et al. (2011), SAT-based probabilistic membership filters Weaver et al. (2014); Schaefer (1978); Douglass et al. (2015); Herr et al. (2017), machine learning applications Hinton (2002); Eslami et al. (2014), or simply estimating the ground-state entropy of a degenerate system. In addition, having multiple solutions to a given problem might allow for the inclusion of constraints in a post-processing step. Here we demonstrate that Monte Carlo methods paired with cluster updates can result in algorithms that asymptotically sample ground-states fairly.

Classical Monte Carlo heuristics based on thermal annealing are known to almost uniformly sample all ground-state and low-lying excited state configurations Moreno et al. (2003); Wang et al. (2015). Studies of three-dimensional diluted Ising antiferromagnets in a field and three-dimensional Ising spin glasses show that parallel tempering Monte Carlo Hukushima and Nemoto (1996) is more efficient than simulated annealing Kirkpatrick et al. (1983) at finding spin-glass ground-state configurations with near-equal probability Moreno et al. (2003); com (a). Isoenergetic cluster moves (ICM) Zhu et al. (2015), related to Houdayer’s cluster updates Houdayer (2001), introduced for Ising spin glasses significantly speed up thermalization on quasi-two-dimensional topologies, such as D-Wave’s Chimera graph. The combination of low-temperature parallel tempering (PT) Monte Carlo and the rejection-free isoenergetic cluster moves (PT+ICM) allow for a wide-spread sampling of search space and help escape local minima separated by large energy barriers. Here we demonstrate that isoenergetic cluster moves paired with parallel tempering Monte Carlo (PT+ICM) enhance the fair sampling of ground-state configurations for spin-glass problems better than the previous PT gold standard. We illustrate the approach using two-dimensional Ising spin glasses on a square lattice, as well as the Chimera graph. Higher-dimensional problems can be embedded in lower-dimensional graphs where PT+ICM is more efficient via, e.g., minor embedding Choi (2008, 2011).

The paper is organized as follows. In Sec. II we introduce a quality metric for fair sampling, as well as a detailed description of a fair-sampling algorithm using ICM. Following that, we present numerical results in Sec. III for both PT, as well as PT+ICM, and introduce an algorithm to approximate the number of degenerate states for highly-degenerate problems. We conclude with a discussion of our results.

II Model and Algorithm

To illustrate the improved sampling of PT+ICM over PT, we start with an Ising spin-glass model on a nonplanar Chimera graph Bunyk et al. (2014). Its nonplanar topology makes finding ground states of random Ising spin glasses worst-case NP-hard. The Hamiltonian for the spin-glass model is given by

[TABLE]

where $s_{i}\in\{\pm 1\}$ are Ising spins and the couplers $J_{ij}$ are drawn for this study from three discrete distributions: $\{\pm 1,\pm 2,\pm 4\}$, $\{\pm 5,\pm 6,\pm 7\}$ and $\{\pm 1\}$). The couplers are selected based on the range of ground-state degeneracy we can handle with our high-performance computing cluster; i.e., the less symmetries between the different coupler values, the smaller the ground-state degeneracy.

II.1 Assessing optimal sampling

Suppose $n$ is the total number of times that ground states are found for an instance with ground-state degeneracy $G$. The probability distribution for finding any particular ground-state configuration follows a binomial distribution. For theoretically perfect sampling, if $p=1/G$ is the probability of finding a state and $q=1-p$ is the probability of failure in a given trial, then the expected number of successes in $n$ trials is $e=np$ and the variance of the binomial distribution is $\sigma^{2}=npq$. Therefore, the theoretical relative standard deviation given by sampling a finite set of random uncorrelated numbers $Q_{\rm th}$ is given by

[TABLE]

Assuming that the states are uncorrelated (which is a safe assumption for large $G$), an algorithm is said to be optimal (sampling fairly) if the numerical relative standard deviation of the frequency of ground-state configurations $Q_{\rm num}$ determined experimentally is close or equal to the theoretical value $\sqrt{(G-1)/n}$ (or $Q_{\rm num}/Q_{\rm th}=1$). In practice, $Q_{\rm num}$ for any algorithm is almost always greater than the theoretical value $Q_{\rm th}$, due to a limited number of measurements via e.g., limited computing resources.

II.2 PT+ICM for fair sampling

Our implementation of PT+ICM for sampling purposes can be summarized as follows:

Run $N_{T}$ replicas of the system at a range of temperatures {$T_{1},T_{2},...,T_{N_{T}}$}, with each set consisting of $M=4$ copies of the system at the same temperature, thus $4\times N_{T}$ copies of the system with the same disorder are randomly initialized. 2. 2.

$N_{\rm sw}$ iterations are performed, each iteration consisting of one Monte Carlo sweep, a parallel tempering update, and an isoenergetic cluster move (for the lowest $N_{\rm hc}$ temperatures). 3. 3.

For the first $N_{\rm sw}/2$ iterations, keep track of the lowest energies for the four replicas at the lowest temperatures. 4. 4.

After $N_{\rm sw}/2$ iterations, the lowest energies $E_{1}$, $E_{2}$, $E_{3}$, and $E_{4}$ for the four replicas with the lowest temperatures are compared, and if $E_{1}=E_{2}=E_{3}=E_{4}$, the ground-state energy has been found with high confidence. Once this is the case, configurations at this energy are recorded, as well as their frequency for the remaining $N_{\rm sw}$/2 updates.

There is no guarantee that any solution obtained by this heuristic method is the true optimum, or that we have found all configurations that minimize the Hamiltonian. However, we choose to make sure each configuration achieves a minimum number of $50$ hits in order to increase our confidence that all accessible ground states have been found. The simulation parameters are shown in Table 1.

III Numerical results

To test whether PT+ICM can sample ground-state configurations with near-equal probabilities, we multiply the numerical relative standard deviation $Q_{\rm num}$ by $\sqrt{n}$ and plot $Q_{\rm num}\sqrt{n}$ as a function of the ground-state degeneracy $G-1$. Note that $Q_{\rm th}\sqrt{n}$ is the square root of the ground-state degeneracy $G-1$, and therefore the function $Q_{\rm th}\sqrt{n}=\sqrt{G-1}$ is a straight line in logarithmic scale for both the horizontal axis ($G-1$) and the vertical axis ($Q_{\rm th}\sqrt{n}$).

Figure 1 shows $Q_{\rm num}\sqrt{n}$ and $Q_{\rm th}\sqrt{n}$ as a function of the ground-state degeneracy $G-1$ for different spin-glass instances on a Chimera graph. As mentioned in the previous paragraph, the quantity $Q_{\rm num}\sqrt{n}$ is almost always greater than $Q_{\rm th}\sqrt{n}$ due to limited computational resources com (b). However, an algorithm samples optimally if the data from the numerical relative standard deviation are close to the theoretical line. It is clear that the data for PT+ICM (blue/darker color) are closer to a straight line than the data for PT (red/lighter color), and the discrepancy between PT+ICM and PT seems to become greater as the system size increases.

In Fig. 2 we plot the median ratio $Q_{\rm num}/Q_{\rm th}$ as a function of the system size $N$ for spin-glass problems on a Chimera lattice. We emphasize that when the ratio becomes unity an algorithm samples optimally. The data show that PT+ICM (blue squares) performs better than PT (red circles) and that the improvement is more significant with increasing system size. In this work the temperature set for the simulation is specifically optimized for $N=1152$. Large median ratios $Q_{\rm num}/Q_{\rm th}$ for smaller system sizes are due to the choice of temperature set. The statistical error bars are determined by a bootstrap analysis using the following procedure: For each system size $N$ and $N_{\rm sa}$ disorder realizations, a randomly selected bootstrap sample of the $N_{\rm sa}$ disorder realizations is generated. The median ratio $Q_{\rm num}/Q_{\rm th}$ is computed with this random sample. We repeat this procedure $N_{\rm boot}=1000$ times for each system size to obtain an average and error bar using these $N_{\rm boot}=1000$ data points.

In addition to studying how fair sampling behaves with increasing system size, we also investigate how the quality of fair sampling is related to ground-state degeneracy and plot $Q_{\rm num}/Q_{\rm th}$ as a function of ground-state degeneracy for $N=800$ variables. Figure 3 suggests that the more ground-state configurations, the easier to sample all ground-state configurations with near-equal probabilities. This is not surprising because a large ground-state manifold makes the algorithm easier to explore the configuration space. We do emphasize, however, that in cases where the ground-state degeneracy is exponentially large and with limited resources only a subset of minimizing configurations is accessible, and for these the sampling improves, the more configurations are present. Furthermore, careful examination of instances with the same system size and ground-state degeneracy suggests that the $Q_{\rm num}/Q_{\rm th}$ ratio is closely related to the Hamming distances between ground-state configurations. It is shown in Fig. 4 that the Sidon instances Katzgraber et al. (2015); Zhu et al. (2016) where $J_{ij}\in\{\pm 5,\pm 6,\pm 7\}$ with large Hamming distances between the ground-state configurations tend to have a high $Q_{\rm num}/Q_{\rm th}$ ratio compared to those with small Hamming distances between the states. Here, PT+ICM achieves more equiprobable sampling with large Hamming distances. Figure 5 shows two examples of ground-state configurations with different Hamming distances on a Chimera graph with $N=128$. PT+ICM’s cluster updates allow nonlocal moves in the energy landscape, therefore reducing $Q_{\rm num}/Q_{\rm th}$ for instances with large Hamming distances between the ground-state configurations.

Figure 6 shows $Q_{\rm num}\sqrt{n}$ and $Q_{\rm th}\sqrt{n}$ as a function of the ground-state degeneracy $G-1$ for different spin-glass instances on a two-dimensional square lattice. Similar to the Chimera graph case, the data using PT+ICM (blue/dark color) are closer to the theoretical optimality line than the data using PT (red/light color), and the discrepancy between PT+ICM and PT becomes larger as the system size increases. In Fig. 7, the median ratio $Q_{\rm num}/Q_{\rm th}$ again demonstrates that PT+ICM is superior to PT in this case for square lattices.

III.1 Estimating the ground-state degeneracy

We also develop an approximate method to count the number of ground-state configurations based on the fair sampling capabilities of PT+ICM and compare the results to exact methods Galluccio et al. (2000) for a handful of configurations. Counting problems Gomes et al. (2008) typically ask how many solutions exist for a given instance and belong to complexity class of #P. This approximate method exploits the fact that if one can sample ground states uniformly then one can obtain a reasonable order-of-magnitude estimate of the ground-state degeneracy. Our renormalization-inspired approach works as follows:

Compute the ground-state energy $E_{0}$ for a fixed number of Monte Carlo sweeps (see above). 2. 2.

Sample the number of ground states $G_{0}$ for the full system for a fixed number of Monte Carlo sweeps. 3. 3.

Iteratively restrict the number of free variables (i.e., those that are not restricted) and estimate the ratio

[TABLE]

for a fixed number of Monte Carlo sweeps. 4. 4.

Repeat until the system size is small enough to be able to compute the number of ground state configurations $G_{\rm final}$ exactly, e.g., via enumeration. 5. 5.

Multiply the product of ratios by the exact count of ground-state configurations to estimate the number of ground states for the full system via

[TABLE]

We compare results of this approximate method to exact counts on a two-dimensional square lattice with bimodal coupling constants $J_{ij}\in\{\pm 1\}$. Simulation parameters and results are shown in Table 2. The renormalization-based estimates agree with the exact ground-state degeneracy within error bars.

IV Conclusions

We have demonstrated that PT+ICM—parallel tempering Monte Carlo with isoenergetic cluster moves—samples ground-state configurations fairly and is an ideal method for applications where a pool of diverse solutions is needed. We also find that degeneracy and Hamming distances between different ground-state configurations are closely related to the relative standard deviation of frequency with which the ground states are found, namely: ground states with large degeneracy and small Hamming distances have a lower relative standard deviation of frequency. It will be interesting to exploit near-uniform sampling for model counting Wei and Selman (2005) and SAT filter construction Weaver et al. (2014); Douglass et al. (2015) in the future.

Acknowledgements.

We would like to thank Ruben S. Andrist, Hamid Khoshbakht, Martin Weigel and Salvatore Mandrà for fruitful discussions. We especially thank Martin Weigel for providing the exact number of ground-state configurations on two-dimensional square lattices using the Vondrak code, and Ruben S. Andrist for rendering Fig. 5. H.G.K. acknowledges support from the National Science Foundation (Grant No. DMR-1151387) and would like to thank Banh Mi for providing the necessary motivation for this research. We thank the Texas Advanced Computing Center (TACC) at The University of Texas at Austin and Texas A&M University for providing HPC resources. Part of 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 purposes notwithstanding any copyright annotation thereon.