Hyperoptimized approximate contraction of tensor networks for rugged-energy-landscape spin glasses on periodic square and cubic lattices

Discrete combinatorial optimization problems and glassy materials with quenched disorder are both often modeled with classical spin Hamiltonians such that the interaction energies between the spins are disordered. The simplest such “spin glass” model is the disordered classical Ising model:

where $i$ and $j$ are graph vertex labels; $J_{ij}$ is a random, real scalar for the interaction energy between spins at vertices $i$ and $j$; and $s_{i}=\pm 1$. Obtaining good solutions to discrete optimization problems and also understanding the low-temperature behavior of glassy systems then entails obtaining the low-energy (a.k.a. near-optimal) and lowest-energy (a.k.a. optimal) configurations of models like Eq. (1) ¹¹1The definition of “near-optimal” is problem-specific. For example, for Max-Cut problems it is in some cases desirable to have solutions within $\sim 12\%$ of optimal and in other cases within $\sim 6\%$ [3]..

Specialized hardware platforms for obtaining near-optimal and optimal configurations of spin glasses are under active development [2, 3, 4, 5, 6]. Although such machines display an average time complexity that is superpolynomial for obtaining optimal solutions [3] (but see Refs. [7, 8, 9, 10] for the theoretical possibility of average polynomial time complexity), the hope is that the smaller constant factors in the scaling that are possible with such machines will one day yield superior performance compared to classical CPU-based methods for industrial optimization problems at scale. A recent milestone of particular note toward this goal is the experimental confirmation by D-Wave, in the context of a 5,000-qubit simulation of a three-dimensional quantum spin glass, of a scaling advantage of quantum annealing over both simulated annealing and (Monte Carlo) simulated quantum annealing for the solution quality (i.e., spin configuration energy distance from the ground state) versus the annealing time [11]. Regarding computing near-optimal solutions (a.k.a. approximate optimization), a recent work [12] showed that quantum annealing can obtain an average polynomial time complexity for solutions within 1% error from optimal of a quasi-two-dimensional classical spin glass, and further that it can do so with a scaling exponent that is smaller than that of the leading classical heuristic (parallel tempering Monte Carlo with isoenergetic cluster moves [13]).

But for both exact and approximate optimization, there is much room left to improve such specialized hardware platforms. For example, densely connected Ising Hamiltonians are the typical case in industrial applications [14], so the sparse connectivity between qubits in the current D-Wave hardware requires an embedding of the dense problem, which incurs a quadratic overhead in the number of spins [14] that dramatically impacts its scaling performance [15]. As another example, the quantum approximate optimization experiments in Ref. [16] show solution quality that remains constant with increasing problem size for spin glass problems with geometry that is native to their superconducting-qubit processor’s planar geometry, but decreasing solution quality with increasing problem size for geometrically non-native problems that are embedded (a.k.a. compiled) into their processor’s geometry. Therefore, further development of classical algorithms for exact and approximate optimization remains relevant, as they can compute larger scale and more complex problems than current specialized machines, and also serve as a way of benchmarking the performance of the specialized machines as they are further developed.

As for classical digital algorithms, many are based on Markov Chain Monte Carlo (MCMC) [17], which performs an annealing-based stochastic search over the low-temperature free-energy landscape in configuration space for the global minima and also local minima that are close in energy to the global minima. However, in various spin-glass cases the ruggedness of the landscape is such that these methods become exponentially slower with increasing system size [18, 19, 20]. There are some proposals to accelerate MCMC methods with assistance from machine learning, but Ref. [19] shows that these proposals fail in the presence of rugged energy landscapes. In fact, Ref. [19] concludes that rugged-energy-landscape spin glass problems may be intrinsically difficult for all machine learning algorithms, and we are not aware of any counterexamples thus far. In fact, the difficulty that MCMC-based methods face in dealing with rugged energy landscapes in the configuration space of spin glass phases has been suggested as indicative of a place to look for quantum advantage [21].

On the other hand, there are classical digital algorithms that are based on simulation of dynamical systems, such as the simulated coherent Ising machine [22, 23] and simulated bifurcation machine [24, 25], but these also amount to a search over the free-energy landscape in configuration space for optimal and near-optimal minima, and can therefore also suffer a steep performance drop (in terms of time to solution) with increasing landscape ruggedness, as observed for example in the supplementary materials of Ref. [23].

For industrial applications it is often necessary to obtain near-optimal solutions to combinatorial optimization problems at very large system sizes ²²2It is often the case in industrial contexts that the optimal solution is not necessary; reduced cost solutions are often sufficient., making polynomial time complexity highly desirable. It is important then to find numerical methods that can approximately solve rugged-energy-landscape classical spin glasses in polynomial time. Tensor networks [27, 28] offer a route to approximately sampling, in polynomial time, the low-temperature Boltzmann distribution of classical spin glasses without traversing the free energy landscape in configuration space, and may thus offer a way to bypass the problems presented by rugged energy landscapes without the use of quantum hardware. This route is now possible due to the following key developments: First, the pioneering work of Nishino [29, 30, 31] showed how to use tensor-network algorithms for homogenous two- and three-dimensional classical lattice models. Then, Murg et al. presented a way in Ref. [32] to use tensor networks to deterministically perform efficient evaluation of partition functions of finite, inhomogeneous two-dimensional classical lattice models with open boundaries (although they did not apply their approach to a spin glass). This latter approach is also reviewed in Ref. [33]. Their idea entails a lattice of low-rank tensors that reflects the geometry of the physical lattice; an approximation of the partition function is obtained in polynomial time by approximately contracting all of the tensors in a row-to-row (or column-to-column) fashion by treating the rows (or columns) as matrix product states (MPSs) or matrix product operators (MPOs), which are one-dimensional tensor networks. The approximation arises due to setting an upper bound on the size of the indices of the intermediate tensors that arise during the contraction. Since the computational complexity of this approximate method is independent of the elements in the tensors, it is independent of the ruggedness of the energy landscape. The subsequent work of Ref. [34] develops the idea of Murg et al. further by demonstrating how to sample actual spin configurations from the low-temperature Boltzmann distribution of spin glasses via such approximate contractions.

We note some limitations of the work in Ref. [34]: It is only able to natively (i.e., without embedding) handle spin glasses on graphs that are either planar or quasi-planar (e.g., the chimera graph), which prevents a native treatment of spin glass problems with non-planar geometry. Though non-planar problems can be embedded into quasi-planar graphs, this entails a polynomial overhead in the number of spins [14]. Further, while Ref. [34] demonstrates good solution quality with small-index-size tensors for large, rugged-landscape spin glasses that are native to the chimera graph, they do not test against problems with non-planar geometry that have been embedded into the chimera geometry. This is significant because tests carried out in Ref. [35] show that non-natively representing the Boltzmann distribution of a complete-graph model by using an MPS requires the size of the tensors of the MPS to grow rapidly with system size if good solution quality is desired. Also, the above-mentioned experiment [16] that found decreasing solution quality with increasing problem size for geometrically non-native problems that are embedded into the native planar geometry of a superconducting-qubit processor makes it doubtful that the quasi-planar tensor-network method in Ref. [34] would yield good solutions in polynomial time for problems that are not natively planar or quasi-planar.

Ref. [36] presents a different tensor network approach for sampling from the partition function, based on tropical algebra, that can natively handle non-planar geometries and that guarantees the optimal solution with a time to solution that is independent of landscape ruggedness. However, the time complexity of this exact approach is (unsurprisingly) superpolynomial in the system size, and due to the use of tropical algebra there is as yet no established way of incorporating approximations into that approach to reduce the time complexity to polynomial. The recent work in Ref. [37] presents a still different tensor network method that is tailored to the exact partition function computation of the p-spin model (which can exhibit rugged energy landscapes) and can also natively handle non-planar geometry, yet again in some regimes has superpolynomial time complexity.

Thus it is important to find ways of extending the polynomial-time approach of Ref. [34] to be able to natively handle non-planar geometries. Strategies for approximately contracting tensor networks with arbitrary geometry do exist [38, 39, 40], and thereby provide a potential route for achieving this. Here we initiate the exploration of this route by testing the ability of the approximate tensor-network contraction method in Ref. [39] to sample the low-temperature Boltzmann distributions of spin glasses on the following non-planar geometries: square and cubic lattices with periodic boundary conditions ³³3Ref. [39] shows that even in cases where an MPS-based contraction, such as the one used in Ref. [34], may be applied, the method of Ref. [39] can significantly outperform MPS-based contraction in terms of both memory and time costs.. We use planted-solution models whose instances span distinct classes of energy landscape ruggedness. For the square lattice we find the solution is always optimal except for in the most rugged class of instances, where we find an energy error of at most 3%. For the three ruggedness classes of the cubic lattice, we find an energy error of at most 12% in the most rugged class, and 0% and about 3% in the other two. These results are all obtained with the same polynomial time complexity ⁴⁴4The (assumed) distinction between P and NP does not preclude the existence of heuristic polynomial-time numerical methods for rugged-energy-landscape spin glasses—methods that have no performance guarantee regarding solution quality but in practice yield good solutions at least in some instances [3].

We follow the general idea in Ref. [34] for representing the Boltzmann distribution $p(\textbf{s})\sim\textrm{exp}[-\beta H(\textbf{s})]$, where s is a spin configuration vector, as a contraction of a network of low-rank tensors; this gives access to any (unconditional or conditional) marginal probability of the system via trivial modifications of the tensor network. Details of our tensor network formulation are given in Appendix A. At a given value of $\beta$, our method computes the marginal probability distribution for a single spin by contracting a modified version of the original tensor network, updates the original tensor network by projecting that spin onto its most probable state, then computes the (conditional) marginal probability distribution of the next spin with the projected version of the original tensor network, and serially iterates this procedure for all remaining spins. In this way the Boltzmann distribution for the entire spin glass at the chosen value of $\beta$ is sampled (though the sampling is not exact, due to the compression described below). We repeat this at a series of values of $\beta$ and choose the lowest energy configuration as the best solution (limited numerical precision precludes the possibility of using arbitrarily large $\beta$).

Ref. [34] performs a more sophisticated sampling via a branch and bound method, but we opt for the simpler strategy above since our focus here is only on providing evidence that hyperoptimized approximate tensor network contraction [39] can, with a fixed polynomial time complexity, readily obtain solutions that are within about 10% energy error from optimal for at least some non-planar spin glass models with rugged energy landscapes.

During the contraction of any tensor network, the size of the indices of the intermediate tensors can grow exponentially. The method in Ref. [39], which we employ in this work, can automatically limit, via compression, the index size to a maximum bond dimension ($\chi$). For short-range spin models on lattices, using a constant value of $\chi$ limits the memory complexity for the tensor contractions to be constant in the number of spins ($N$). The leading order memory complexity of the entire algorithm is therefore only that of the initial tensor network itself, which is $\mathcal{O}(N)$ for short-range models (see Appendix A for an example). For a fixed $\chi$, the time complexity in terms of $N$ is set by the number of tensor contractions. For short-range spin models on lattices, there are $\mathcal{O}(N)$ tensor contractions for each single-spin marginal. Computing a spin configuration for a given $\beta$ with the method described above for all $N$ spins therefore has a time complexity $\mathcal{O}(N^{2})$. We note that the computations at different values of $\beta$ and $\chi$ are independent, and can therefore be done in parallel.

As explained in Ref. [34], while increasing $\beta$ theoretically reduces thermal fluctuations and therefore the $\chi$ that is needed to achieve a given solution quality (larger $\chi$ yields less compression and therefore more accuracy in principle), limited numerical precision means that in practice there is a finite value of $\beta$ that is optimal at any given value of $\chi$ (to see this, consider that the numerical range of the Boltzmann weights in Eq. (2) in Appendix A monotonically increases with increasing $\beta$). Finite numerical precision also means that increasing $\chi$ does not monotonically decrease the error at a fixed $\beta$. Thus the values of $\chi$ and $\beta$ that optimize the solution quality are not known a priori. But as mentioned above, our purpose here is only to provide evidence that the general scheme is sufficiently performant in the context of rugged, non-planar spin glasses to be of interest for combinatorial optimization. Here we say that this means that solutions within about $10\%$ of optimal should appear even at small values of fixed $\chi$, since computational cost monotonically increases with increasing $\chi$. Hence we usually set $\chi=4$ and test a few values of $\beta$ to show that multiple solutions with energy error $\lesssim 10\%$ do appear.

While the actual time cost of any tensor network contraction can depend strongly on the order of the tensor contractions, the contraction algorithm in Ref. [39] that we employ automatically chooses an optimal or near-optimal contraction order. This search for the contraction order is done only once for a given lattice geometry, and not for every spin glass instance, so we do not include its cost in the computational complexity of the algorithm. We implement the algorithm with the Python libraries quimb [43] and cotengra [44].

Our testing uses the planted-solution Ising-model instances generated by the Chook library [45] with the tile planting method [46] on the square lattice and cubic lattice, both with periodic boundaries in all directions; in both cases the Hamiltonian takes the form in Eq. (1). Since the solutions are planted, the groundstate energy is exactly known.

The simulations are done one instance at a time on an Apple M2 Ultra processor with 24 cores (16 performance cores and efficiency cores) and 128 GB of RAM.

We have investigated the ability of approximate tensor network contraction with small maximum bond dimension to yield optimal or near-optimal (i.e., $\lesssim 10\%$ energy error) solutions in polynomial time to at least some types of non-planar, short-range spin glass problems with rugged free energy landscapes. For planted-solution Ising spin glasses generated with tile planting, with a constant maximum bond dimension of $\chi=4$, and therefore fixed polynomial time and memory complexity, we found no substantial reduction in solution quality with increasing system size between $400$ and $2304$ spins in the toroidal-square-lattice case, and between $64$ and $216$ spins in the toroidal-cubic-lattice case. For the most rugged subclass of instances in the toroidal-square-lattice case, the average best solution quality (over multiple instances) remained at $\sim 1\%$ with increasing system size, while it remained at $\sim 10\%$ for the most rugged subclass in the toroidal-cubic-lattice case. As far as solution quality, these results are encouraging given that the sampling method that we used was quite crude; it is likely that a more sophisticated sampling procedure, such as the one in Ref. [34], would yield a substantial reduction in the minimum energy error for some of the problem instances without increasing $\chi$. At the same time, we emphasize that the method here does not have a strict guarantee on solution quality (i.e., it is categorized as a heuristic method). Regarding system size, though the time complexity of the method is quadratic, we were unable to reach beyond $6\times 6\times 6$ for the cubic-lattice due to the actual simulation time required; further development of the method is needed in order to be able to use tensor-network representations of partition functions to approximately optimize the much larger cubic-lattice sizes addressable with certain other methods.

For instance, the recent D-Wave quantum annealing experiment in Ref. [11] also yielded an approximate optimization of a cubic-lattice spin-glass model, but at a much larger size of over 5,000 spins while achieving an error of about 2$\%$ on a much shorter time scale of microseconds. A direct comparison between their results and the ones here is not possible since they did not assess the computational time complexity (they used a fixed system size), nor did they test the performance of their method against different levels of energy-landscape ruggedness, which in the literature has become synonymous with computational hardness [47, 19]; it is known that if a method successfully performs approximate optimization for some spin glass models at very large sizes (e.g., 8,000 spins on a cubic lattice [48]), it may still perform poorly on energy landscapes that are more rugged [19]. But we can make a rough comparison between our results and those of Ref. [11] as follows: According to Ref. [46], the most rugged subclass ($F_{6}$) for the cubic lattice model that we test is orders of magnitude more computationally difficult than the random $\pm J$ instances that are tested in Ref. [11]. Considering this, our obtained error of $\sim$10% for subclass $F_{6}$ (Figs. 5 (e-f)) suggests that our method (or a further development thereof) may perform competitively in terms of solution quality (i.e., relative error) for the spin glass instances of Ref. [11]. This is further suggested in Figs. 5 (a-d), which show that when $L=6$ our method achieves an error of about 3% for the $F_{22}$ subclass and an error of 0% in the $F_{42}$ subclass. Also, Fig. 6 (b) reveals an average error of $\sim 5\%$ in the $F_{6}$ subclass when $L=4$ and $\chi=32$. But as mentioned above, it is unlikely that the method here could approximately optimize a cubic lattice as large as 5,000 spins in a practical amount of time, and in that sense one can say that the D-Wave results are superior. At the same time, the method here has a quadratic time complexity, and further development of this method may therefore eventually yield something competitive regarding reachable system sizes; the fact that small bond dimension ($\chi$) yields good results suggests that the physics is highly local, and ongoing followup work [49] suggests that this might be leveraged to dramatically improve simulation time and thereby reach much larger system sizes. On the other hand, there is theoretical hope that the time complexity of quantum annealing may become universally polynomial [50, 8, 9, 10], which could help quantum annealing reach still larger system sizes. Also noteworthy is that Ref. [51] demonstrates that cubic spin glasses in the form of minor embeddings that are larger than what can fit on a single quantum annealing processor can also be approximately optimized with quantum annealing by a hybrid quantum-classical approach that outperforms classical simulated annealing at a fixed system size. Thus, the situation regarding competition between methods for optimization of cubic-lattice Ising spin glasses remains fluid.

It is also worthwhile to juxtapose our square-lattice results with those of the D-Wave quantum annealing experiment in Ref. [12], which uses a quasi-two-dimensional model whose interaction graph is equivalent to a honeycomb lattice with additional non-planar, short-range bonds. Our results (Fig. 3) show 0% error except in the most rugged subclass ($C_{2}$), where we achieve an average error of $\lesssim$1%, and our time complexity is always quadratic. The error-corrected quantum-annealing and parallel-tempering Monte-Carlo data in Ref. [12] are for error thresholds at and very close to 1%, and always have sub-quadratic but super-linear time complexity. The Sidon-set [52] spin glass instances that they used in Ref. [12] are considerably computationally harder than $\pm J$ spin glass instances, but we do not know which of the subclasses of the present square-lattice model it is closest to in terms of hardness.

Regarding the possibility of addressing more complex classical spin glasses than quasi-two-dimensional Ising and cubic-lattice Ising, the tensor-network approach investigated here is amenable to beyond-binary variables and beyond-quadratic interactions, and can easily incorporate local fields. Further, Ref. [53] provides a constructive proof that square- and cubic-lattice classical Ising spin glasses faithfully capture the low-energy configuration space of any classical spin system with only a polynomial overhead in the number of spins (see also Ref. [54] regarding mapping between K-SAT and the cubic-lattice classical Ising spin glass). When dealing with such square- and cubic-lattice spin glasses that encode more general spin glasses whose interaction graphs are highly non-local, it is likely that the method investigated here would suffer from a significant bias in the sampling due to the finite bond dimension (the non-local nature of the original model that is encoded into the lattice model would likely require very large bond dimension for the present method to sample the lattice model approximately fairly). However, some work has demonstrated that unbiased sampling of Boltzmann distributions of classical spin lattices can be achieved with tensor-network contraction by combining it with Markov Chain Monte Carlo [55, 56, 57, 58]. In particular, Ref. [58] demonstrated this to work well with small bond dimension for a classical square-lattice spin glass at criticality, where the correlation length is very large; it may be that the method would work well at small bond dimension also for planar or cubic-lattice spin glass models that encode spin glasses with highly non-local interaction graphs. It is unclear, however, whether or not the time complexity would be polynomial.

Finally, short-range spin glass models are better models of real materials with quenched disorder than complete-graph models, and therefore are of relevance for condensed matter physics. A prominent example is the classical Edwards-Anderson Ising model on the cubic lattice, the nature of the low-temperature phase of which has remained unresolved due to the rugged-energy-landscape problem [59]. Since the results in this work demonstrate that tensor-network contraction can approximately optimize rugged-energy-landscape cubic spin glasses of moderate size in polynomial time, we hope that the developments here and in related work may eventually contribute toward resolving such fundamental questions.

We acknowledge discussions with Garnet Chan and Yu Tong. AAG also acknowledges correspondence with Sam Reifenstein, Timothée Leleu, Wangwei Lan, and Marek Rams.