Enumeration algorithms for combinatorial problems using Ising machines: When should we stop exploring energy landscapes?

TL;DR

The paper tackles enumerating all optimal or constraint-satisfying solutions in combinatorial problems by leveraging Ising machines as solution samplers. It introduces stopping criteria rooted in coupon-collector theory and formalizes energy-ordered fair and cost-ordered fair samplers, providing theoretical bounds that ensure the probability of incomplete enumeration is at most $\epsilon$. Two algorithms are developed: Algorithm 1 for constraint satisfaction and Algorithm 2 for combinatorial optimization, each with provable guarantees and practical stopping rules; the latter maintains provisional optima and uses a subroutine to enumerate feasible solutions at a given cost. A numerical demonstration on maximum clique enumeration with simulated annealing shows competitive computation times against a branch-and-bound baseline on dense graphs, though sampler fairness can vary in practice; results indicate significant potential for fast, enumeration-enabled decision-making in real-world problems. The work highlights both the benefits and limitations of sampling-based enumeration and points to future directions, such as using alternative samplers with guaranteed fairness and integrating Boltzmann-like sampling to broaden applicability.

Abstract

Combinatorial problems such as combinatorial optimization and constraint satisfaction problems arise in decision-making across various fields of science and technology. In real-world applications, when multiple optimal or constraint-satisfying solutions exist, enumerating all these solutions is often desirable, as it provides flexibility in decision-making. However, combinatorial problems and their enumeration versions pose significant computational challenges due to combinatorial explosion. To address these challenges, we propose enumeration algorithms for combinatorial optimization and constraint satisfaction problems using Ising machines. Ising machines are specialized devices designed to efficiently solve combinatorial problems by exploring the energy landscape of an Ising model. Ising machines typically sample lower-energy solutions with higher probability. Our enumeration algorithms repeatedly perform such sampling to collect all desirable solutions. The crux of the proposed algorithms lies in their stopping criteria for sampling-based energy landscape exploration, which are derived from probability theory. In particular, the proposed algorithms have theoretical guarantees that the failure probability of enumeration is bounded above by a user-specified value, provided that lower-cost solutions are sampled more frequently and equal-cost solutions are sampled with equal probability. Many physics-based Ising machines are expected to (approximately) satisfy these conditions. As a demonstration, we applied our algorithm using simulated annealing to maximum clique enumeration on random graphs. We found that our algorithm enumerates all maximum cliques in large, dense graphs faster than a conventional branch-and-bound algorithm specifically designed for maximum clique enumeration. These findings underscore the effectiveness and potential of our proposed approach.

Figures (8)

• Figure 1: Pseudocode of Algorithm 1. The function SAMPLE is a fair sampler of feasible solutions. The definition of $\kappa_1(\epsilon)$, which appears in line 6, is provided in Eq. \ref{['eq:kappa1-def']}. The failure probability of Algorithm 1 is theoretically guaranteed to be less than the user-specified failure tolerance $\epsilon$ (see Theorem \ref{['theorem:algorithm-1']} in Appendix \ref{['appx:algorithm-1']}).
• Figure 2: An illustration of a sampling process in Algorithm 1. Each circle represents a sample generated by an Ising machine, with different colors indicating different solutions. The value above each circle indicates the energy of the corresponding solution. In this example, the ground state energy is 0.0, and thus light-blue and green circles represent feasible solutions (i.e., desirable solutions). In contrast, circles of other colors correspond to infeasible solutions, which are discarded during the sampling process, as indicated by "$\times$" marks. The numbers below the feasible solutions indicate the sample count $\tau$. The deadlines for $m = 2$ and $3$ were calculated as $\lceil m\ln(m\kappa_1/\epsilon) \rceil$ with $\epsilon = 0.01$.
• Figure 3: Pseudocode of Algorithm 2. The function SAMPLE is a cost-ordered fair sampler of feasible solutions. The definition of $\kappa_2(\epsilon)$, which appears in line 5, is provided in Eq. \ref{['eq:kappa2-def']}. The failure probability of Algorithm 2 is theoretically guaranteed to be less than the user-specified failure tolerance $\epsilon$ (see Theorem \ref{['theorem:algorithm-2']} in Appendix \ref{['appx:algorithm-2']}).
• Figure 4: An illustration of a sampling process in Algorithm 2. Each circle represents a sample generated by an Ising machine, with different colors indicating different solutions. The "$\times$" marks indicate the rejection of the samples. The value above each circle indicates the energy of the corresponding solution. Here, for simplicity, all samples are assumed to be feasible, and the energy of each feasible solution is set to its cost value in the original problem. In this example, the ground state energy is 0.0, and thus light-blue and green circles represent the optimal solutions. The numbers below the accepted solutions indicate the sample count $\tau$. Note that the sample count is reset when the current minimum cost value $\theta$ is updated. The deadlines for $m = 2$ and $3$ were calculated as $\lceil m\ln(m\kappa_2/\epsilon) \rceil$ with $\epsilon = 0.01$. The sample sequence is the same as that in Fig. \ref{['fig:algorithm-1-example']}.
• Figure 5: Computation times required to enumerate all maximum cliques in random graphs with different graph densities $D$ and different numbers of vertices. The blue circles indicate the computation times of Algorithm 2 using SA. For each data point, the mean computation time of all successful cases out of 100 independent runs was calculated. For graphs with density $D=0.75$ and more than 355 vertices [the right-hand side of the vertical dotted line in panel (c)], where the conventional exact algorithm did not terminate even after 10 days and thus the true maximum cliques are unknown, we treat runs in which all largest cliques found across the 100 runs were obtained as "successful" cases for the purpose of estimating the computation time. The relative standard errors had a mean of 0.01 and a maximum of 0.07; thus, the error bars (not shown) are shorter than the diameter of the blue circles. The sky-blue dashed lines represent linear fits to the computation times indicated by the blue circles for graphs with more than 250 vertices. The red diamonds indicate the computation times of the conventional algorithm (Bron--Kerbosch combined with the enhancements by Tomita--Tanaka--Takahashi and Carraghan--Pardalos). For cases with computation times shorter than seven days, the average of 10 independent runs was taken. For cases with computation times longer than seven days, only one run was conducted to evaluate the computation time. The relative standard errors had a mean of 0.01 and a maximum of 0.33; thus, the error bars (not shown) are shorter than the size of the red diamonds. The light-pink dashed lines represent linear fits to the computation times indicated by the red diamonds for graphs with more than 250 vertices.

Theorems & Definitions (14)

• Lemma 1
• proof
• Lemma 2
• proof
• Lemma 3
• proof
• Lemma 4
• proof
• Theorem 1
• proof