Solving Spin Glasses with Optimized Trees of Clustered Spins

TL;DR

This paper introduces an algorithm that constructs optimized trees of clustered spins for efficient optimization and thermal sampling of spin glasses, outperforming existing methods on various problem classes.

Contribution

The paper proposes a novel algorithm combining tree construction and thermal sampling for spin glasses, improving efficiency and scalability over prior approaches.

Findings

Outperforms existing algorithms on benchmark problems.

Effectively balances tree size and sampling complexity.

Demonstrates versatility across different graph connectivities.

Abstract

We present an algorithm for the optimization and thermal equilibration of spin glasses - or more generally, cost functions of the Ising form $H=\sum_{\langle i j\rangle} J_{ij} s_i s_j + \sum_i h_i s_i$, defined on graphs with arbitrary connectivity. The algorithm consists of two repeated steps: i) the optimized construction of a random tree of spin clusters on the input problem graph, and ii) the thermal sampling of the generated tree. The randomly generated trees are constructed so as to optimize a balance between the size of the tree and the complexity required to draw Boltzmann samples from it. We benchmark the algorithm on several classes of problems and demonstrate its advantages over existing approaches.

Figures

Equations

$$H=\sum _{⟨ij⟩}{J}_{ij}{s}_i{s}_j+\sum _i{h}_i{s}_i$$

$$E_{\text{inter}}({\mathbf{s}}_i,{\mathbf{s}}_j)=\sum _{k,k^{\prime }}{J}_{i(k),j(k^{\prime })}{s}_{i(k)}{s}_{j(k^{\prime })},$$

$$E_{\text{self}}({\mathbf{s}}_i)=\sum _k{h}_{i(k)}{s}_{i(k)}+\sum _{k<k^{\prime }}{J}_{i(k),i(k^{\prime })}{s}_{i(k)}{s}_{i(k^{\prime })}.$$

$$E(t_{{\mathbf{s}}_i})$$

$$W(t_{{\mathbf{s}}_i})={\text{e}}^{-\beta E_{\text{self}}({\mathbf{s}}_i)}\prod _j\left[{\text{e}}^{-\beta E_{\text{inter}}({\mathbf{s}}_i,{\mathbf{s}}_j)}W(t_{{\mathbf{s}}_j})\right],$$

$$P(t_{{\mathbf{s}}_i})=\frac{W(t_{{\mathbf{s}}_i})}{{\sum }_{t_{{\mathbf{s}}_i}}W(t_{{\mathbf{s}}_i})}.$$

$$P(t_{{\mathbf{s}}_i})=P({\mathbf{s}}_i,t_{/{\mathbf{s}}_i})=P(t_{/{\mathbf{s}}_i}|{\mathbf{s}}_i)P({\mathbf{s}}_i).$$

$$P({\mathbf{s}}_i)=\frac{P({\mathbf{s}}_i,t_{/{\mathbf{s}}_i})}{P(t_{/{\mathbf{s}}_i}|{\mathbf{s}}_i)}=\frac{{\sum }_{t_{/{\mathbf{s}}_i}}W(t_{{\mathbf{s}}_i})}{{\sum }_{t_{{\mathbf{s}}_i}}W(t_{{\mathbf{s}}_i})},$$

$$P({\mathbf{s}}_i)=\frac{W({\mathbf{s}}_i)}{{\sum }_{{\mathbf{s}}_i^{\prime }}W({\mathbf{s}}_i^{\prime })},$$

$$W({\mathbf{s}}_i)$$

$$C_{\mathbf{t}}=\sum _{⟨{\mathbf{s}}_i,{\mathbf{s}}_j⟩}{2}^{|{\mathbf{s}}_i|}\times 2^{|{\mathbf{s}}_j|}\approx 2^{2⟨|\mathbf{s}|⟩}(|\mathbf{t}|-1)$$

$$⟨|\mathbf{s}|⟩=\frac{1}{|\mathbf{t}|}\sum _i|{\mathbf{s}}_i|$$

$$F=\frac{C_{\mathbf{t}}}{2^{{\sum }_i|{\mathbf{s}}_i|}}=\frac{2^{2⟨|\mathbf{s}|⟩}(|\mathbf{t}|-1)}{2^{|\mathbf{t}|⟨|\mathbf{s}|⟩}},$$

Full text

I Introduction

Many problems of theoretical and practical relevance consist of searching for the global minimum of intricate cost functions defined over discrete configuration spaces. These so-called combinatorial optimization problems are not only notoriously hard to solve, but are also ubiquitous, appearing in a wide range of diverse fields such as machine learning, materials design, software verification, flight-traffic control, constraint satisfaction problems and logistics, to name a few diverse examples Wolsey and Nemhauser (1999); Papadimitriou and Steiglitz (2013). It is no surprise then that the design of fast and practical algorithms to solve these problems is an important challenge in many areas of science and technology.

A class of problems that is known to be particularly challenging to solve is that of spin-glasses Young (1998)—disordered, frustrated spin systems, so intricate that specialized hardware has been built to simulate them Belletti et al. (2008, 2009); Baity-Jesi et al. (2014). The spin glass, or Ising model, cost function is usually written as

[TABLE]

defined over a set of binary variables, the Ising spins $\{s_{i}=\pm 1\}$. The spins sit at the vertices of an interaction graph defined by the sets of given parameters $\{J_{ij}\}$ and $\{h_{i}\}$. Here, $\sum_{\langle ij\rangle}$ denotes a sum over all the edges of the connectivity graph of the problem. Finding the minimal cost and spin assignments that produce it, or using physics terminology, the ground state energy and corresponding spin configurations for problem instances belonging to the above model, is generally known to be a combinatorially intractable (or, NP-hard) optimization problem Barahona (1982). Furthermore, a significant portion of real-life hard discrete optimization problems can be cast as spin glasses Lucas (2014), making the development of novel practical algorithms with significantly shorter runtimes and/or superior scaling with problem size the holy grail of optimization.

The problem of finding the lowest costs of spin-glass Hamiltonians, Eq. (1), can be and is often generalized to an analogous problem in statistical mechanics, where one is asked to draw spin configurations from the Boltzmann distribution $\mathcal{B}(T)$ of that cost function for a given temperature $T$—or expressed differently, to sample the configuration space of the problem such that the probability of drawing any given spin configuration ${\mathbf{s}}=(s_{1},\ldots,s_{N})$ is proportional to its Boltzmann weight $W_{\mathbf{s}}={\text{e}}^{-\beta E_{\mathbf{s}}}$ where $\beta=1/T$ and $E_{\mathbf{s}}$ is the cost associated with the configuration. In the zero temperature limit $T\to 0$, Boltzmann sampling reduces to the standard optimization variant of the problem wherein only ground state assignments have nonzero weights.

The Boltzmann sampling of spin glasses at nonzero temperatures has numerous immediate practical applications in of itself in various areas (e.g., in machine learning Ackley et al. (1985)). Furthermore, and perhaps more importantly, Boltzmann sampling is known to serve as a useful subroutine in optimization algorithms where only the absolute lowest cost is sought. The observation that the thermalization, or equilibration, of intricate cost functions at sequentially decreasing temperatures is beneficial for finding optimal assignments was first made by Kirkpatrick et al. Kirkpatrick et al. (1983) who, based on this observation, pioneered the now-celebrated approach commonly referred to as ‘simulated annealing’ (SA). Today, the methods of choice for studying general spin-glass problems are variations of SA. One scheme is ‘population annealing’ Hukushima and Iba (2003) in which simulated annealing is combined with Boltzmann weighted differential reproduction within a population of replicas to sample equilibrium states. Another technique, that we describe in more detail in Sec. II.3, is ‘parallel tempering’ (PT) Hukushima and Nemoto (1996); Marinari (1998) in which multiple copies of the problem are equilibrated in parallel at different temperatures and spin configurations at adjacent temperatures are regularly swapped.

At the most basic level, the equilibration of spin glasses at any given temperature $T$ consists of the repeated updating of a stored spin configuration ${\mathbf{s}}$ where usually one calculates the energy change $\Delta E=E_{{\mathbf{s}}_{\text{new}}}-E_{\mathbf{s}}$ associated with the flipping of a single spin while keeping the values of all other spins fixed and accepting the move with a probability that satisfies a detailed balance, namely, $P_{{\mathbf{s}}\to{\mathbf{s}}_{\text{new}}}/P_{{\mathbf{s}}_{\text{new}}\to{\mathbf{s}}}=W_{{\mathbf{s}}_{\text{new}}}/W_{\mathbf{s}}$. Here, ${\mathbf{s}}$ denotes the current spin configuration and ${\mathbf{s}}_{\text{new}}$ is the proposed new configuration. The above condition is sufficient to ensure the eventual proper thermalization of the system. For spin glasses, updates consisting of single spin-flips (SSFs) are often found to be grossly inefficient. This is mainly because the system often finds itself trapped at local minima (or metastable states) of the cost function, usually surrounded by high energy barriers. The probability for a spin configuration to ‘jump over’ such barriers is exponentially suppressed in the barrier height $\Delta E$. On the other hand, moves consisting of the thermal sampling of larger subsystems, i.e., multiple spin-flips, which correspond to potentially higher jumps, are generally very costly to implement, as the complexity of calculating the energy changes and Boltzmann weights associated with multiple spin-flips grows exponentially with the size of the spin cluster that is to be flipped.

Certain subgraph structures, however, can nonetheless be sampled, or optimized, more efficiently. The best example is trees—subgraphs that do not contain any cycles. Trees are known to give rise to equilibration techniques, such as belief propagation Yedidia et al. (2003), which allow for a more efficient Boltzmann sampling, or ground state optimization, than that of single spin updates Selby (2014); Hamze and de Freitas (2004); Decelle and Krzakala (2014). The computational complexity associated with the Boltzmann sampling of trees of spins scales only linearly with their size. However, optimization of spin glasses via trees of single spins (TSSs) induced on input problem connectivity graphs is very often found to be inefficient compared to single spin-flip techniques. In part, this is because arbitrarily connected spin glass graphs may contain many short cycles, and thus severely restrict the size of the induced trees. Moreover, the computational overhead associated with the generation of random trees on the instance connectivity graph plays a role in diminishing the scaling advantage of using them Decelle and Krzakala (2014).

In this work we propose a novel algorithm for the thermalization and optimization of spin glasses. The algorithm aims to adequately address the pitfalls of existing approaches discussed above—namely, the inefficiency of single spin-flips and the impracticality that is often associated with the sampling of trees of single spins. The algorithm we propose here is based on heat-bath updates performed on optimized trees of spin clusters induced on the input problem graph. The algorithm consists of two basic repeated steps. In the first step, a tree of clusters of spins is induced on the spin glass graph in a manner which optimizes a certain balance between the size of the subgraph that is to be thermalized and the complexity of doing so. As a next step, a configuration from the Boltzmann distribution $\mathcal{B}(T)$ of the tree is drawn (while keeping the spins not belonging to the tree fixed). As we demonstrate, a methodical grouping together of spins into clusters generally allows for the construction of much larger trees on the input graph—a process which in turn allows for the efficient thermalization of large subgraphs, thereby generally speeding up the entire equilibration process.

This paper is organized as follows. In Sec. II we discuss the two steps of the algorithm, namely, the processes for generation and subsequent sampling of trees of spin clusters. We also illustrate how one can use the complexity of sampling a tree of clustered spins to construct a figure of merit with which optimal subgraph structures can be generated on the input problem graph and then thermalized. In Sec. III we present the results of several benchmarking tests comparing the performance of the proposed algorithm, as an optimizer as well as for thermal equilibration, to existing techniques—specifically, single spin-flip and random single spin tree methods. Section IV is devoted to a discussion of the results and the applicability of the method to optimization problems of practical relevance.

II Thermal annealing with optimized trees of spin clusters

For reasons that will become clear shortly, we begin our description of the algorithm with the second of its two steps, namely, the subroutine for the Boltzmann sampling of a tree of clustered spins. The complexity of this subroutine, which we calculate here, will play a part in the design of the first step of the algorithm, within which the tree structures that are to be thermalized in the second step are constructed.

II.1 Boltzmann sampling of a tree of clustered spins

Given an input problem of the Ising type, Eq. (1), we define a cluster as a set of spins of the input problem that may or may not be coupled among themselves. A tree of spin clusters (or a TOSC) is a subgraph of the input problem graph that has the structure of a tree—that is, each node on the tree is a cluster of spins, and edges between neighboring nodes exist if and only if Ising interactions between their respective spins exist. Figure 1(a) illustrates a tree of clustered spins.

We now discuss an algorithm for drawing a spin configuration on a TOSC with Boltzmann probabilities of a given inverse temperature $\beta$. Fixing the values of all the input-problem spins not belonging to the tree, interactions with outside spins will henceforth be treated as fixed external fields. Figure 1(b) illustrates a tree of clustered spins embedded in a larger input problem graph. For what follows, we shall denote tree nodes (i.e., spin clusters) or configurations thereof by $\mathbf{s}_{i}=\{s_{i(1)},s_{i(2)},\ldots\}$ where $i$ labels the node, and $s_{i(k)}$ denotes the individual spins comprising the node.

For brevity, we shall denote the combined interaction energy between spins in neighboring nodes $i$ and $j$ as

[TABLE]

and the energy associated with the interaction of spins within a node, combined with the energy associated with external fields (including the interaction energy with fixed spins not in the tree), as

[TABLE]

We consider now a tree of clustered spins with a root node $\mathbf{s}_{i}$. Let us denote the set of spins comprising the tree by $\mathbf{t}_{\mathbf{s}_{i}}$ and configurations of spins on the tree by $t_{\mathbf{s}_{i}}$. Denoting the set of children nodes of $\mathbf{s}_{i}$ by $c(\mathbf{s}_{i})$, and $\mathbf{s}_{j}\in c(\mathbf{s}_{i})$ as the children nodes, we have $\mathbf{t}_{\mathbf{s}_{i}}=\{\mathbf{s}_{i}\}\cup_{j}\mathbf{t}_{\mathbf{s}_{j}}$; i.e., the set of clustered spins forming a tree is the union of the root node with all the (sub-)tree sets whose root nodes are $\mathbf{s}_{j}$. The energy $E(t_{\mathbf{s}_{i}})$ of any given tree configuration $t_{\mathbf{s}_{i}}$ can thus be written as

[TABLE]

From the above definition, it follows that the weight associated with a configuration $t_{\mathbf{s}_{i}}$ can be written as

[TABLE]

where the goal of the algorithm is to generate a spin configuration $t_{\mathbf{s}_{i}}$ with a probability that is proportional to $W(t_{\mathbf{s}_{i}})$; that is,

[TABLE]

To generate a thermal configuration over the tree $\mathbf{t}_{\mathbf{s}_{i}}$ with a Boltzmann probability, we will first calculate the probability $P(\mathbf{s}_{i})$ of assigning a specific value to $\mathbf{s}_{i}$ (equivalently, specific values to the spins of $\mathbf{s}_{i}$). After $\mathbf{s}_{i}$ is fixed, the problem immediately decouples to the thermalization of a set of disconnected (sub-)trees $\mathbf{t}_{\mathbf{s}_{j}}$ for which the process will be repeated.

To calculate the probability $P(\mathbf{s}_{i})$ of assigning a specific value $\mathbf{s}_{i}$ to the root node, we shall use the relation $P(A\cap B)=P(A|B)P(B)$, from which it follows that the probability for drawing a configuration of spins on the tree can be written as the product

[TABLE]

where $t_{\not{\mathbf{s}}_{i}}$ is any configuration of spins over the set $\mathbf{t}_{\mathbf{s}_{i}}$ with the value of $\mathbf{s}_{i}$ fixed. Rearranging, we find that:

[TABLE]

where $\sum_{t_{\mathbf{s}_{i}}}$ denotes a sum over all spin configurations on the set $\mathbf{t}_{\mathbf{s}_{i}}$, and $\sum_{t_{\not{\mathbf{s}}_{i}}}$ denotes a sum over all configurations over the set $t_{\mathbf{s}_{i}}$ with the value of the root node set to $\mathbf{s}_{i}$. Rewriting $\sum_{t_{\mathbf{s}_{i}}}=\sum_{{\mathbf{s}_{i}}}\sum_{t_{\not{\mathbf{s}}_{i}}}$, the above becomes:

[TABLE]

where we have defined the weight of a node $\mathbf{s}_{i}$ as:

[TABLE]

We thus find that the probability of assigning a specific value to any clustered spin $\mathbf{s}_{i}$ can be calculated from its weight, which is obtained from the weight of the tree of which $\mathbf{s}_{i}$ is the root. After fixing the root node $\mathbf{s}_{i}$, the probabilities of its children $\mathbf{s}_{j}$ can similarly be determined and their values subsequently fixed.

To assign a configuration to all nodes in a tree, we follow these two steps: i) Advancing from the tree leaves (that is, from nodes with no children) to root, we calculate the weights of all the possible node configurations. The calculation of the weight of a node, Eq. (5), is determined by its self-energy, the interaction energy with its children and their weights. ii) Having obtained the weights of the various nodes, we begin assigning configurations to the tree nodes, starting from the root node advancing toward the tree leaves. Assigning a configuration to a node becomes possible once the configuration of its parent node has been fixed, at which point the node becomes a root node for a subtree, and the node’s self energy is amended to include its interaction with its (now fixed) parent node.

Since the evaluation of the weight of every node configuration requires the calculation of the node’s self energy, its interaction with its children and the weights of its children, the complexity, or number of elementary steps, required for assigning values to the spins of a tree of clustered spins, is linear in the number of edges (or nodes) in the tree. Specifically, it requires $2^{|\mathbf{s}_{i}|}\times 2^{|\mathbf{s}_{j}|}$ evaluations (where $|\mathbf{s}_{i}|$ is the number of spins in $\mathbf{s}_{i}$), for every pair of neighboring nodes $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$:

[TABLE]

where

[TABLE]

is the average number of spins per node and $|\mathbf{t}|$ is the size of the tree; i.e., the number of tree nodes (here, $|\mathbf{t}|-1$ is the number of edges).111In the above, we keep implicit the linear dependence of the complexity of the calculation on the degree (or number of neighbors) of the spins in the tree.

II.2 Locally optimal random trees of spin clusters

Having computed the complexity of drawing a configuration from the Boltzmann distribution of an Ising model defined on a tree of clustered spins (or in the $T\to 0$ limit, the complexity of finding an optimal configuration on that tree), we turn to addressing the question of generating optimal TOSCs for a given input problem—TOSCs from which it would be most advantageous in terms of overall equilibration time to draw Boltzmann configurations. On the one hand, the larger the subgraphs induced on the input problem graph are, the lower the number of ‘tree flips’ required to thermalize the input graph. On the other hand, the complexity associated with thermalizing trees with large clusters may become prohibitively costly. To capture the needed balance between tree size and the complexity of drawing a thermal Boltzmann sample from it, we define the simple figure of merit $F$ — the ratio of the complexity of drawing a thermal configuration from a tree $C_{\mathbf{t}}$ to the number of configurations on it:

[TABLE]

which captures the relative ease with which one can draw a thermal configuration on the tree out of $2^{|\mathbf{t}|\langle|\mathbf{s}|\rangle}$ possible configurations. The smaller $F$ is, the more merit there is to its thermalization.

We are now in a position to outline the first step of the TOSC algorithm in which a random TOSC is generated ‘on-the-fly’ on the input graph, minimizing (locally) the figure of merit $F$. This first step is based on sequentially adding spins to a subgraph, re-evaluating $F$ at every stage, until a local minimum of $F$ is found, at which point the process is ended.

The generation of a random TOSC on a given spin glass connectivity graph is carried out as follows. At first, a random spin that will function as a root node for the TOSC is selected. As a next step, its neighboring spins (or spins of distance one from the root) are added to the root, one by one, each forming a new neighboring tree node. The formation of a new node trivially lowers the value of $F$ and is thus always desired. If a spin that is to be added to the constructed tree is found to close a cycle, i.e., if it is connected to two or more different nodes on the tree, adding it as a new node as is, is prohibited. At this point, two alternatives are considered: i) The new spin is added to the TOSC by ‘contracting’ the cycles formed due to its inclusion. This is done by unifying nodes in a manner that eliminates the cycles. Figure 2 provides illustrations of cycle contraction, as the algorithm prescribes. ii) In the second alternative, the spin is not added to the TOSC being constructed and is taken out of the pool of addable spins. To determine which alternative is more beneficial, we evaluate $F$ for each scenario and pick the alternative with the smaller value. The process continues in a similar manner with the addition of spins that are neighbors of the neighbors of the root (equivalently, spins of distance two from the root) and so forth, until there are no additional spins to add. Figure 3 illustrates the formation of a random TOSC embedded in a given input problem.

The full algorithm thus consists of two steps: the first in which a random TOSC is generated by optimizing the benefit of thermalizing large subgraphs, and the second which actually draws a thermal sample from the tree.

II.3 Parallel tempering with TOSC

A way to further speed up the TOSC algorithm is to embed it in a parallel tempering (PT) scheme. In PT, one considers $N_{T}$ independent copies of the system running in parallel at different temperatures, $T_{1}<T_{2}<\ldots<T_{N_{T}}$.222Here, we choose a temperature grid of the PT simulations consisting of $N_{T}=30$ temperatures. Temperatures with indices $i=13,14,\ldots,N_{T}$ were evenly distributed in the range $0.21\leq T_{i}\leq T_{\text{max}}=1.632$, while lower temperatures in the range $T_{\mathrm{min}}=0.045\leq T_{i}\leq 0.2$ (indices $i=1,2,\ldots,12$) V. Martin-Mayor and I. Hen (2015). Copies with neighboring temperatures regularly attempt to swap their temperatures with probabilities that satisfy detailed balance Sokal (1997). In this way, each copy performs a temperature random-walk. At high temperatures, free-energy barriers are easily overcome, allowing for a global exploration of configuration space. At lower temperatures on the other hand, the local minima are explored in more detail. We will consider a variant of the algorithm where a single TOSC is generated for all replicas at the beginning of each ‘temperature sweep,’ and during the sweep random samples are drawn on the tree for each copy independently, based on its temperature.

III Algorithm benchmarking

We now present the results of several benchmarking tests performed on the proposed algorithm. We compare the typical runtimes of three parallel tempering-based optimization algorithms: (i) single spin-flip Metropolis algorithm (SSF), (ii) parallel tempering with randomly generated trees of single spins (TSS), and (iii) the TOSC algorithm.

Here, the typical runtimes for any set of parameters is defined as the median time to reach a minimizing configuration (i.e., a ground state), or thermal equilibrium, over 100 randomly generated instances with the specified parameters. All the algorithms were run on one core of a $3.5$ GHz 6-Core Intel Xeon E5 processor.

III.1 Optimization of trees of fixed-size spin clusters

The first class of random Ising problems we test is one where the underlying connectivities of the input problems themselves are randomly generated trees of spin clusters of fixed-size $C$, wherein all the spins within a cluster are connected by edges, and edges also exist between all spin pairs that belong to neighboring clusters. An example of a random $20$-node input problem with $C=4$ spins per cluster (and $80$ spins in total) is given in Fig. 4. The coupling strengths $J_{ij}$ on each edge are chosen to be $\pm 1$ with equal probability, and the external fields $h_{i}$ are set to zero.

For this type of problem, an a priori knowledge of the structure of the input problem allows for the devising of an optimization algorithm that takes full advantage of the underlying tree structure of the problem. In this case, the runtime scales linearly with the number of clusters for any fixed $C$. In Fig. 5, we show the typical runtimes of the SSF, TSS and TOSC algorithms as a function of problem size for different values of cluster size $C$. As is immediately evident, the typical runtimes of all algorithms scale polynomially with problem size. However, for both SSF and TSS the power (the slope) grows rapidly with increasing cluster size $C$, whereas the TOSC runtime scales linearly with problem size irrespective of $C$, similar to the ideal algorithm. For $600$-spin problems, this amounts to speed-ups of up to four orders of magnitude of TOSC as compared to the other algorithms.

It is also instructive to look at the average ratio of tree size $\langle|\mathbf{t}|\rangle$ to problem size $N$ of the trees generated by TSS and TOSC. This is shown in Fig. 6. While the TSS tree covers less and less of the input problem with growing cluster size [Fig. 6(a)], the TOSC trees tend to cover the entire input graph in the limits of large problem sizes [Fig. 6(b)].

III.2 Thermalization of trees of fixed-size spin clusters

To test the performance of the TOSC algorithm as a thermalizer rather than an optimizer, that is, to ascertain its ability to draw configurations from the Boltzmann distributions of given Hamiltonians at given temperatures, we next compare the typical runtimes of SSF, TSS and TOSC to equilibrate trees of fixed-size spin clusters. To that aim, we utilize the concept of PT mixing time $\tau$ Fernandez et al. (2009); Alvarez Baños et al. (2010); Fernandez et al. (2013); V. Martin-Mayor and I. Hen (2015); Marshall et al. (2016) using it as the figure of merit for thermalization. The mixing time of a PT simulation may be thought of as the average time it takes a PT replica to fully traverse the temperature mesh, indicating equilibration of the simulation. In PT, copies with neighboring temperatures regularly attempt to swap their temperatures with probabilities that satisfy detailed balance Sokal (1997). In this way, each copy performs a temperature random-walk. At high temperatures, free-energy barriers are easily overcome, allowing for a global exploration of configuration space. At lower temperatures on the other hand, the local minima are explored in more detail. A ‘healthy’ PT simulation requires an unimpeded temperature flow and so the total length of the simulation should be longer than the temperature mixing time Fernandez et al. (2009); Alvarez Baños et al. (2010).

In Fig. 7 we plot as a function of problem size the typical ‘time to thermalization,’ which we take here to be the time it takes all PT replicas to traverse the temperature grid from the lowest temperature to the highest and back four times. Similar to the optimization results of the previous section, we find that typical thermalization runtimes scale polynomially with problem size for all three algorithms. However, while for SSF and TSS the power (the slope) seems to grow with increasing cluster size $C$, the TOSC runtime scales much more mildly, indicating a strong scaling advantage (the errors in the reported powers are of the order of the second decimal place).

III.3 Spin glasses on Chimera graphs

The problem classes we consider next are random spin glasses whose underlying connectivity graphs are of the Chimera type: two-dimensional arrays of unit cells of eight spins with a $K_{4,4}$ bipartite connectivity Choi (2008, 2011). A 2-cell by 3-cell subgraph of a Chimera lattice is shown in Fig. 8. While the Chimera graph is two-dimensional in nature, it is also non-planar and as such gives rise to difficult spin-glass problems Barahona (1982). Chimera graphs have become the focus of much attention in recent years in the context of quantum annealing-based optimization due to the commercial availability of prototypical quantum annealing optimizers of spin glasses whose quantum bits are coupled with a Chimera connectivity Johnson et al. (2011); Berkley et al. (2013); Bunyk et al. (2014).

Here, we measure the typical runtimes to reach a ground state configuration as recorded by the SSF, TSS and TOSC algorithms on spin glass instances with Chimera connectivities where the couplings are chosen randomly from either $J_{ij}\in\{\pm 1\}$ (‘Range $1$’ instances) or $J_{ij}\in\{\pm 1,\pm 2,\pm 3\}$ (‘Range $3$’). The results are summarized in Fig. 9. The first and second panels show the speed-up, exhibited by the milder slope, gained by TOSC as compared to SSF and TSS for both Range $1$ (a) and Range $3$ (b) classes (the errors in the reported slopes are of the order of the second decimal place). Figure 9(c) depicts the average coverage of of the trees of TSS (about $50\%$) and TOSC (approximately $75\%$). Also shown is the TOSC average cluster size $\langle|\mathbf{s}|\rangle$.

III.4 Algorithmic scaling with problem hardness

We next study the manner in which the performance of the TOSC algorithm correlates with ‘instance hardness.’ Since instances with large $\tau$ are harder to equilibrate, we follow Refs. Fernandez et al. (2009); Alvarez Baños et al. (2010); Fernandez et al. (2013); V. Martin-Mayor and I. Hen (2015); Marshall et al. (2016), and use as a measure for the hardness of a random spin glass instance, the mixing time $\tau$ of a single spin-flip PT simulation.

To benchmark the TOSC algorithm, we generate about $10^{6}$ random instances on an $N=512$-spin Chimera graph and measure the mixing time of each instance V. Martin-Mayor and I. Hen (2015); Marshall et al. (2016). As a next step, we group together instances with similar classical hardness, i.e., similar mixing times, $10^{k}\leq\tau\leq 3\cdot 10^{k}$ for $k=3,4,5,6$ and $7$. For each such ‘generation’ of $\tau$, we randomly pick $100$ representative instances for the benchmarking of the algorithm (only 14 instances with $k=7$ were found). The results are presented in Fig. 10 which shows on a log-log scale the typical runtimes of the SSF, TSS and TOSC algorithms as a function of mixing time. As is evident, the susceptibility of both SSF and TSS to classical hardness scales approximately linearly with mixing time with a slope slightly below $1$. On the other hand, the TOSC algorithm is found to be much less sensitive to instance hardness, with a slope that is about three times smaller, indicating a clear qualitative advantage over the other algorithms.

IV Summary and conclusions

We developed a novel algorithm for the optimization and sampling from the Boltzmann distribution of Ising spin glasses of arbitrary sizes and connectivities. The algorithm is designed to take advantage of the connectivity graph of the spin glass by thermalizing randomly generated subgraphs of the input problem—specifically, trees of spin clusters (TOSCs) constructed to optimally balance between the size of the thermalized subgraph and the complexity of doing so.

Benchmarking the TOSC algorithm against single spin-flip PT and single spin random tree PT, we showed that the TOSC algorithm provides qualitative scaling advantages on all problem classes that have been tested, both in terms of scaling with problem size as well as scaling with instance hardness. In light of the scaling advantages of the algorithm discussed here and the generality of its scope, we believe that the algorithm will prove to be a useful optimization tool for potentially many classes of problems of practical relevance (once these are cast as spin glasses).

It would be of interest to explore additional and possibly more efficient algorithms for generating random optimal subgraphs based on the figure of merit introduced above, or alternate figures of merit, with which optimal subgraphs are induced and then thermalized on the input problem graph. It is also worth noting that for problem classes that share a single, specific underlying connectivity graph, the generation of globally optimal trees of clustered spins may also be useful, as these may be generated in advance and used in a manner similar to the random structures suggested here. Furthermore, the TOSC algorithm may also be combined with other efficient approaches in lieu of PT such as the Houdayer algorithm Houdayer (2002).

Acknowledgements.

We thank Victor Martin-Mayor for many insightful discussions on the topics addressed in this paper, and Tameem Albash for useful comments and suggestions.