A Continuous Energy Ising Machine Leveraging Difference-of-Convex Programming

TL;DR

The paper addresses the challenge of solving large-scale Ising ground-state problems by relaxing binary spins to continuous variables and introducing an attraction term that biases solutions toward binary configurations. It formulates the resulting Hamiltonian as a difference of convex functions and develops two DCP-based solvers, DOCH and ADOCH, with provable convergence properties and low per-iteration cost suitable for GPU acceleration. Empirical results show that DOCH/ADOCH outperform state-of-the-art solvers across small to ultra-large problem sizes, including a 10^7-spin fully connected model with tens of trillions of couplings solved in hours on GPUs. The approach provides a scalable, cooling-free alternative for large-scale combinatorial optimization with broad implications for high-performance Ising-based computation.

Abstract

Many combinatorial optimization problems can be reformulated as finding the ground state of the Ising model. Existing Ising solvers are mostly inspired by simulated annealing. Although annealing techniques offer scalability, they lack convergence guarantees and are sensitive to the cooling schedule. We propose solving the Ising problem by relaxing the binary spins to continuous variables and introducing an attraction potential that steers the solution toward binary spin configurations. A key property of this potential is that its combination with the Ising energy produces a Hamiltonian that can be written as a difference of convex polynomials. This enables us to design efficient iterative algorithms that require a single matrix-vector multiplication per iteration and provide convergence guarantees. We implement our Ising solver on a wide range of GPU platforms, from edge devices to high-performance computing clusters, and demonstrate that it consistently outperforms existing solvers across problem sizes ranging from small ($10^3$ spins) to ultra-large ($10^8$ spins).

Figures (26)

• Figure 1: Overview of our Ising solver for a $2$-spin system. We relax the $2^2 = 4$ discrete energy values $(\mathcal{E}_0, \mathcal{E}_1, \mathcal{E}_2, \mathcal{E}_3)$ to a smooth energy landscape $\mathcal{E}(x_1,x_2)$ defined over the continuous space $(x_1,x_2) \in \mathbb{R}^2$. This energy is combined with an attraction potential $\mathcal{A}(x_1,x_2)$ to form the Hamiltonian $\mathcal{H} = \mathcal{E} + \mathcal{A}$. The local minima of this Hamiltonian, denoted by green dots $(\mathcal{H}_0, \mathcal{H}_1, \mathcal{H}_2, \mathcal{H}_3)$, serve as candidate ground states. With suitable parameter tuning, the minimizer $(x^*_1, x^*_2)$ of the Hamiltonian produces a binary vector $\mathrm{sign}(x^*_1, x^*_2)$ that closely matches the true ground state $\boldsymbol{s}^* = (s^*_1, s^*_2)$ of the original Ising energy $\mathcal{E}$.
• Figure 2: Shape of the attractor.(a) Plot of the one-variable attractor for different $\alpha, \beta$ values. The global minimizers at $\pm (\alpha/\beta)^{1/2}$ are shown. (b) Plot of the attractor in two variables for $\alpha = \beta = 1$, showing the four global minimizers at $(1, 1), (1, -1), (-1, 1)$ and $(-1, -1)$.
• Figure 3: Optimization trajectory for a $2$-spin system. We consider the anti-ferromagnetic model $J_{12} = J_{21} = -1$, whose ground states are $(1,-1)$ and $(-1,1)$. Top row: Surface plots of the attractor $\mathcal{A}$, Ising energy $\mathcal{E}$, and the Hamiltonian $\mathcal{H}=\mathcal{A}+\mathcal{E}$. Bottom row: Corresponding contour plots for $\alpha = 1$ and $\beta = 2$. In the bottom-right plot, the white curves represent the optimization trajectories of our Ising solver, starting from various initial points (indicated by green dots). Depending on the initialization, the iterates converge to one of the two degenerate ground states (red dots) after five iterations.
• Figure 4: Evolution of Ising energy with runtime for the $10^3$-spin SK model.(a) Comparison of our solvers (DOCH and ADOCH) with state-of-the-art Ising solvers: Simulated Annealing (SA) isakov2015inagaki2016, ballistic Bifurcation Machine (bSB) goto2021, Simulated Coherent Ising Machine (SimCIM) goto2021, and Spring Ising Algorithm (SIA) jiang2024. Solid curves show the best energy values achieved across $100$ independent runs (Supplementary Notes 6 and 7 for the parameter settings). The dashed black line indicates the lowest energy obtained using the GW-SDP goto2019inagaki2016. (b), (c) Histograms of $100$ Ising energies, each obtained after $1000$ iterations and with different initializations. Results obtained on a NVIDIA RTX 3050 GPU laptop.
• Figure 5: Ising energy versus runtime on the $K_{2000}$ benchmark.(a) Comparison of our solvers (DOCH and ADOCH) with state-of-the-art Ising solvers. Solid curves show the mean Ising energy across $100$ independent runs (Supplementary Notes 6 and 7 for the parameter settings). The dashed black line indicates the lowest energy obtained using the GW-SDP goto2019inagaki2016. (b), (c) Histograms of $100$ Ising energies, each obtained after $1000$ iterations and with different initializations. Results obtained on a NVIDIA RTX 3050 GPU laptop.

Theorems & Definitions (26)

• Proposition 1
• proof
• Proposition 2
• proof
• Proposition 3
• proof
• Proposition 4
• proof
• Theorem 1
• proof