Neural Ising Machines via Unrolling and Zeroth-Order Training

TL;DR

The paper addresses NP-hard Ising/Max-Cut optimization by learning a compact, data-driven update rule for iterative Ising-machine dynamics using a two-layer MLP whose parameters evolve in time via a Fourier-based basis. Training uses a zeroth-order evolutionary optimizer (DAS) to overcome unstable gradients from long recurrent trajectories, yielding NPIM with emergent momentum and annealing-like behavior, while maintaining a small parameter footprint. Across standard benchmarks, NPIM and especially its discrete variant dNPIM achieve competitive or state-of-the-art performance in objective value and time-to-solution on Max-Cut, MIS, and Max-Clique problems, including large G-set graphs up to $N=2000$. The approach offers a scalable, interpretable data-driven heuristic that can adapt to a range of combinatorial optimization tasks, with potential extensions to broader problem classes and hybrid gradient-based training.

Abstract

We propose a data-driven heuristic for NP-hard Ising and Max-Cut optimization that learns the update rule of an iterative dynamical system. The method learns a shared, node-wise update rule that maps local interaction fields to spin updates, parameterized by a compact multilayer perceptron with a small number of parameters. Training is performed using a zeroth-order optimizer, since backpropagation through long, recurrent Ising-machine dynamics leads to unstable and poorly informative gradients. We call this approach a neural network parameterized Ising machine (NPIM). Despite its low parameter count, the learned dynamics recover effective algorithmic structure, including momentum-like behavior and time-varying schedules, enabling efficient search in highly non-convex energy landscapes. Across standard Ising and neural combinatorial optimization benchmarks, NPIM achieves competitive solution quality and time-to-solution relative to recent learning-based methods and strong classical Ising-machine heuristics.

Figures (7)

• Figure 1: a) Diagram showing the flow of information in a neural network Ising machine. The coupling fields (purple), calculated by aggregating the influence of the $N-1$ other variables, are saved from previous iterations. They are then fed into the neural network model and used to decide the new spin variable which is then used to compute the next coupling field. See Secs. \ref{['sec: IM']} and \ref{['sec: IM_MLP']} for a concrete mathematical description. b) High level overview of our method relative to other approaches to CO, inspired by Fig. 15 of Monga2020AlgorithmUnrolling. c) Cartoon depiction of the zeroth-order evolutionary optimization algorithm we use based off of Reifenstein2024DynamicAnisotropicSmoothing. A distribution of parameters (represented by red circles) is evolved over many iterations to close in on an optimal parameter configuration (blue dot).
• Figure 2: Example of training single layer neural network Ising machine. Upper left: the average reward (success rate) of the network with respect to training epoch. The reward starts negative because of an initial bootstrapping phase in which bad trajectories are penalized. Two snapshots of the network parameters are taken and shown in the right two figures: network A at epoch 19, and network B at epoch 99. Lower left: residual Ising energy (difference with best-known solution) is shown as a function of iteration step for both networks. Darker colored trajectories indicate the ground state was found. Bottom middle and right: network weights of network A and network B respectively. Blue and red connections depict negative and positive network weights respectively. Top middle and right: trajectory of $x_i(t)$ variables for network A and network B respectively. Each color represents a different variable of the Ising problem.
• Figure 3: a) Performance (Time to solution) of cNPIM on Sherrington-Kirkpatrick (SK) problem instances. Colored traces show performance with and without fine-tuning showing limited (but nonzero) ability to generalize over problem size. Dotted trace shows baseline Ising machine performance (Chaotic amplitude control Leleu2019Leleu2021Leleu2025). b, e) Scatter plot showing the TTS of 100 random SK problem instances of problem size $N=800$ against that of CAC for cNPIM and dNPIM respectively. c) Success rate for different architectures of cNPIM and dNPIM as a function of total parameter count. The same data is shown in Tab. \ref{['tab: arch_table']}. d) TTS is shown as a function of hardness parameter for the Wishart planted ensemble (WPE) problem instances Hamze2020WPE. Colored traces show cNPIM fine-tuned on different hardness parameters while dotted line shows Ising machine baseline (CAC).
• Figure 4: Hyperparameter sweeps of NPIM on SK instances. a) Final success rate as a function of history length $T_c$ with $D=3$ and $M=3$. b) Success rate as a function of hidden neurons $D$ with $T_c=8$ and $M=3$. c) Success rate as a function of Fourier modes $M$ with $T_c=8$ and $D=3$. Curves compare the continuous (cNPIM) and discrete (dNPIM) variants trained on $N=100$ SK instances. Batch size $B=20$, and the success-rate reward.
• Figure 5: Comparison of temporal basis functions in NPIM. Final success rate for Fourier, Legendre, and Chebyshev bases as a function of degrees of freedom per parameter $M$, with $T_c=8$ and $D=3$ fixed. Left: cNPIM. Right: dNPIM. Networks are trained on $N=100$ SK instances using the success-rate reward, with $400$ epochs, $R=400$ trajectories per epoch, and batch size $B=20$. Performance improves as $M$ increases, and differences between bases are small once sufficient degrees of freedom are available.