Scalable Continuous-Energy Ising Solvers Leveraging Difference-of-Convex Programming

Optimizing the Ising model to its ground state is a central NP-hard combinatorial problem arising across physics, materials science, machine learning, and operations research. Practical solvers often rely on heuristic annealing methods, which can scale effectively but require cooling schedules and parameter tuning and typically provide limited convergence guarantees. At the other end of the spectrum, continuous relaxations such as semidefinite programming offer rigorous guarantees but become computationally prohibitive as problem size grows. We propose a continuous relaxation in which the discrete spin space $-1,1$${}^n$ is extended to ${\mathbb{R}}^n$. A separable quartic attractor replaces the hard constraint; its minimizers are exactly the vertices of a scaled hypercube $-\lambda ,\lambda$${}^n$. After solving the resulting continuous problem, a discrete solution is obtained via a componentwise sign mapping. The sum of the Ising energy and the attractor admits a difference-of-convex polynomial structure. This yields the Difference-of-Convex Hamiltonian (DOCH), which requires one matrix-vector multiplication and a componentwise cube-root update. Using the analytic nature of the Hamiltonian and the Łojasiewicz gradient inequality, we prove that DOCH converges to a critical point in this parameter regime and, under a nonzero-coordinate and Jacobian-stability condition, that the limit is a strict local minimum, enjoying linear convergence. For accelerated DOCH (ADOCH), we prove finite-window envelope descent and an $\mathcal{O}(N^{-1/2})$ best-iterate stationarity bound in the strict DC regime, together with full-sequence convergence when the look-back window is set to zero ($q=0$) from some finite iteration onward. We implement the proposed solvers on GPU platforms ranging from edge devices to high-performance computing clusters. Across MAX-CUT and spin-glass benchmarks, matched-budget comparisons with annealing, physics-inspired, and recent sampling-based solvers show favorable runtime-quality trade-offs from $10^3$ to $10^8$ spins.