Graph Coloring Approach to Solving Sudoku with Oscillatory Neural Networks Filip Sabo and Aida Todri-Sanial
arXiv:2607.15814v1 [cs.LG] 17 Jul 2026
NanoComputing Research Lab, Electrical Engineering Department Eindhoven University of Technology Eindhoven, the Netherlands [email protected], [email protected]
Abstract—Oscillatory Neural Networks (ONNs) present an attractive physics-based computing paradigm rooted in the dynamics of a network of typically fully coupled oscillators aiming to minimize an underlying energy function. In this paper, we propose an ONN-based solver for one well-known constrained combinatorial optimization problem, namely a Sudoku, by formulating the problem as a Graph Coloring problem. By modifying the already existing Graph Coloring solver to a computationally cheaper version and introducing an additional term ensuring the fulfillment of the Sudoku constraints, our solver is shown to significantly outperform the existing HNNand ONN solvers in terms of accuracy. In particular, we are able to achieve nearly flawless accuracies on 4 × 4 as well as rather high accuracies on 9 × 9 Sudoku puzzles for different numbers of unknown digits. Index Terms—Oscillatory Neural Networks, Coupled Oscillators, Graph Coloring, Sudoku.
I. I NTRODUCTION Constrained combinatorial optimization problems (COPs) are ubiquitous in the industry, ranging from scheduling through chip design to the finance sector [1–4]. In particular, NP-hard problems present a special subclass of COPs where the time to find the solution becomes exponential as the problem size increases [5]. Additionally, the current computing paradigm is not suitable for dealing with such problems, as it suffers from the memory–compute bottleneck resulting in high power consumption [6–8]. In this regard, there has recently been a push for alternative computing paradigms. In this paper, we investigate the performance of a physical computing paradigm based on a network of coupled oscillators, called Oscillatory Neural Networks (ONNs) [9]. In contrast to conventional Artificial Neural Networks [10– 12], Oscillatory Neural Networks present an emerging physicsbased computing paradigm rooted in the dynamics of a network of typically fully coupled oscillators aiming to minimize an underlying energy function [9, 13, 14]. In such networks, the information is encoded in the phases of the oscillators. Thanks to their close ties to Hopfield Neural Networks [15, 16] and the Ising model [17, 18], ONNs are well suited for autoassociative memory tasks [19–27] and combinatorial optimizaThis work has received funding from the European Union’s Horizon Europe research and innovation programme, PHASTRAC project under grant agreement No 101092096 as well as European Research Council ERC THERMODON project under grant agreement No. 101125031.
tion problems [14, 28–35]. In this paper, we propose an ONNbased solver for one well-known constrained combinatorial optimization problem, namely a Sudoku. Sudokus embody the renowned digit placement puzzle. In particular, given a grid of size N × N and some known digits, the aim lies in filling out the remaining empty cells of the grid with integers from the range [1, N ] such that no repetition of digits occurs in any row, column, or box of the grid. Thus, Sudokus can be viewed as a constrained combinatorial optimization problem. In particular, Graph Coloring, the problem of assigning each node in a graph a “color” such that a minimal number of colors is utilized and no two adjacent nodes share the same color, offers a well-suited fit for Sudokus. In this paper, we build an ONN-solver tailored to Sudoku’s constraints based on the Graph Coloring problem. Recently, the authors in [28] have developed an ONN-based Sudoku solver by embedding the given Sudoku constraints into the weight matrix. Their method has been shown to improve upon the existing HNN-based Sudoku solver for 4×4, 9 × 9 as well as 16 × 16 puzzles in terms of accuracy [36]. In contrast, in this paper, we present an alternative ONNSudoku solver based on Graph Coloring with an additional term ensuring the fulfillment of the constraints. Our model is found to significantly outperform the existing HNN- and ONN-solver in terms of accuracy on 4 × 4 and 9 × 9 Sudoku puzzles by reaching almost perfect accuracies on 4 × 4 and predominantly high accuracies on 9 × 9 Sudokus for various numbers of unknown digits. This paper is structured as follows. In Section II, we introduce the key concepts behind physics-based computing with Oscillatory Neural Networks (ONNs). Furthermore, Section III outlines our proposed ONN-Sudoku solver based on Graph Coloring together with the benchmarking set-up. Most importantly, our benchmarks on 4 × 4 and 9 × 9 Sudokus in terms of accuracy against the established HNN- [36] and ONN-Sudoku solvers [28] are presented in Section IV. Last but certainly not least, the results are discussed in Section V and the paper is concluded with Section VI. II. BACKGROUND A. Ising model This paper revolves around Oscillatory Neural Networks (ONNs), an emerging physics-based paradigm rooted in the
dynamics of coupled oscillators [9]. In contrast to Artificial Neural Networks [10–12], ONN computing boils down to finding global minima in an underlying energy function [9, 13, 14]. The ONNs are inherently linked to the Ising model [17, 18], a mathematical model of ferromagnetism describing a spin configuration of a system [17]. In particular, the dynamics of the model can be described with the following Hamiltonian: X H=− Jij σi σj . (1) (i,j)
Given a coupling weight matrix Jij , the goal lies in identifying the natural state of the system by minimizing the outlined Hamiltonian through the set of binary spins σi . B. Oscillatory Neural Networks (ONNs) One can bridge the gap between the Ising model and ONNs by mapping the spins σi ∈ {1, −1} to the oscillators’ phases ϕi ∈ {0, π} resulting in the corresponding ONN-Hamiltonian. Moreover, the ordinary differential equation (ODE) describing the dynamics of N oscillators in an ONN can thus be determined simply by following the negative gradient of the ONN-Hamiltonian [14, 19]: 1 X Wij sin (ϕj − ϕi ), (2) ϕ̇i = N j where Wij stands for the coupling between the i-th and j-th oscillator. The ODE in Eq. 2 is often called the Kuramoto model [37]. As already mentioned, the phases (0, π) encode the Booleans (1, −1). However, when running an ONN, it rarely happens that the phases settle at the pre-established binary values 0 and π, but rather at continuous ones. Hence, to force the oscillators into the binary phases, one introduces an additional term scaled by a factor KS into the ODE called Subharmonic Injection Locking (SHIL) which injects the second harmonics into the oscillators [9, 14, 38]: 1 X ϕ̇i = Wij sin (ϕj − ϕi ) − KS sin (2ϕi ). (3) N j III. M ETHODS In this Section, we introduce our approach for solving Sudoku puzzles with Oscillatory Neural Networks. Since the rules of Sudoku imply non-repeating values across all columns, rows, and boxes, Graph Coloring presents a suitable approach for solving this problem. At the end of this Section, we outline the benchmarking process. A. Proposed ONN-Sudoku solver 1) Graph Coloring with ONNs: Given a graph, the aim of the Graph Coloring problem boils down to assigning each node in a graph a color such that no two adjacent nodes share the same color and a minimal number of colors is utilized. Being a part of the 21 Karp NP-complete problems [5], Graph Coloring is a well-known hard combinatorial optimization problem. As shown by the author in [18], all NP-hard problems can be mapped to a corresponding Ising Hamiltonian with Graph Coloring not being an exception.
Furthermore, authors in [35] demonstrated the mapping for some non-binary hard optimization problems to Oscillatory Ising Machines (“Oscillatory Neural Networks”). In particular, Graph Coloring is defined as a Max-K-Cut with the minimal K such that no two adjacent nodes share the same set. In the case of a Sudoku, the number of sets or K is well-known and does not need to be minimized, as it corresponds to the size of the Sudoku. Thus, the desired phases ϕideal for a size K Sudoku read as follows: 2π(k − 1) k ∈ {1, 2, ..., K} . (4) ϕideal = K Hence, for a 9 × 9 Sudoku (K = 9), we get the following mapping of the numbers to the phases: 1→0 2π 2→ 9 4π 3→ 9 2π 4→ 3 8π 5→ (5) 9 10π 6→ 9 4π 7→ 3 14π 8→ 9 16π 9→ 9 Let us now turn our attention towards Max-K-Cut. Given a graph, the Max-Cut problem entails finding a partition between the nodes into 2 sets such that the sum of the cut edges is maximized. In the case of Max-K-Cut, a partition into K sets with the maximum sum of the cut edges shall be found. The authors in [35] solve Max-K-Cut with Oscillatory Neural Networks by extending ODE from Eq. 3 to the following form: 1 X ϕ̇i = − Wij sin (ϕi − ϕj + f (ϕi − ϕj )) N j (6) − KS sin (Kϕi ), with f (∆ϕij ) ensuring that the ideal phases from Eq. 4 produce a vanishing gradient when acted upon with a sine function. Furthermore, the function f reads as follows: K−1 X 2kπ f (∆ϕij ) = (2k − 1)π − K k=1 ( ) 2 (∆ϕij − 2kπ K ) (7) exp − 2σ 2 ( ) 2 (∆ϕij + 2kπ K ) − exp − , 2σ 2 with σ being a tuneable parameter. Fig. 1 depicts the function f (∆ϕij ) from Eq. 7 and the term sin (∆ϕij + f (∆ϕij )) for
2
1 0 Phase difference
1
ij( )
2
ij))
0.5
ij + f(
1.0
3 2 1 0 1 2 3
0.0
sin(
f(
ij)( )
K = 4 and a relatively large σ (= 0.15) to showcase the form of the function.
0.5 1.0
2
1 0 Phase difference
1
ij( )
2
Fig. 1: Function f (∆ϕij ) from Eq. 7 (left) and term sin (∆ϕij + f (∆ϕij )) (right) for K = 4. Clearly, utilizing the model from Eq. 6 has multiple drawbacks. First of all, the function in Eq. 7 has a tuneable parameter σ, which controls the width of the peaks. Second, computing the terms in Eq. 7 is computationally expensive, particularly when done so for multiple thousands of time steps and large problems. However, most importantly, the process of reaching the roots of the term sin (∆ϕij + f (∆ϕij )) is not identical, since the gradients can vary significantly depending on the chosen route. This can lead to certain roots being preferred or in turn avoided completely. Thus, we need an easily computable model without any tuneable parameters and where each root is identical. All of these conditions are satisfied with the following model: K 1 X Wij sin (ϕi − ϕj ) − KS sin (Kϕi ). (8) ϕ̇i = − N j 2 In particular, the first term ensures that the oscillators are 2πn K apart, while the second term guarantees that the oscillators settle at 2πn K values. We utilize Eq. 8 instead of Eq. 6 in our Sudoku approach. Based on Eq. 8, we define an order parameter κ to track how well the oscillators converged to the desired phases as outlined in Eq. 4: iKϕi 1 X . (9) κ= 1 − Im exp N i 2 The order parameter is bounded to the range [0, 1], where values close to 1 indicate that the phases ϕi converged well to the pre-established phases ϕideal , whereas values near 0 hint at a poor convergence to ϕideal . Let us now discuss the form of the weight matrix Wij . 2) Weight matrix: Let us now turn our attention to the weight matrix. First of all, we impose a zero diagonal. Second, thinking of the rules of Sudoku, there should not be any repeating values in each row, column and box. Thus, each cell should only be connected to all cells in its corresponding row, column and box. Since the Graph Coloring problem is defined as a Max-K-Cut with the minimal K such that no two adjacent nodes share the same set, the value connecting any two cells in a common row, column or box should be negative. Without loss of generality, we opt for a value of −1.
Third and most importantly, not all weights should have equal strength, as the knowns should have a stronger influence on the dynamics than the unknowns. This implies that a weight between two unknowns (let us call it wu ) shall be negative and have a smaller magnitude than the weight connecting a known (−1 in our case). Putting all of these conditions together, we arrive at the following form of the weight matrix Wij : −1, i ̸= j, i ∧ j ∈ row, column or box, i ∨ j ∈ known, Wij = wu , i ̸= j, i ∧ j ∈ row, column or box, i ∧ j ∈ / known, 0, else, (10) with the following two conditions: • wu < 0, • |wu | < 1. While the weight matrix in Eq. 10 guarantees that the knowns will have a stronger influence on the unknowns than the unknowns, the ODE does not ensure that the unknowns will have no influence on the knowns. Hence, to produce a vanishing gradient for knowns, we modify the ODE presented in Eq. 8 to the following form: ϕ̇i = − funknown (ϕi )
1 X Wij sin N j
K (ϕi − ϕj ) 2
(11)
− funknown (ϕi )KS sin (Kϕi ), with ( funknown (ϕi ) =
1, ϕi ∈ unknown, 0, else.
(12)
3) Ensuring different phases: While the model presented in Eq. 11 accordingly describes the ONN-ODE for Max-K-Cut, it has a serious shortcoming in regards to Sudoku. Figure 2 showcases an instance of a simple 4×4 Sudoku puzzle with a violation of the rules. 1
2
1
1 2
3 2
2 3
Fig. 2: Example of a 4×4 Sudoku grid with a violation of the rules, where the thick numbers indicate the known values. Clearly, the digit 1 in the upper left hand side corner violates two rules simultaneously, as there is already a known digit 1 in the same column as well as box. However, despite breaking the rules, the ODE presented in Eq. 11 would produce a vanishing gradient, for the digits 1 and 1 are indeed 2πn K apart (n = 0 is also permitted) and both values are one of the valid roots.
In summary, despite the structure, the ODE presented in Eq. 11 allows the ONN to break the Sudoku rules. In order to circumvent this issue, we need an additional term in the ODE. This new term shall kick the oscillators out of a stable state if a violation of the rules is detected. The ODE of the modified proposed approach to solve Sudoku with MaxK-Cut reads as follows: 1 X ϕ̇i = − funknown (ϕi ) Wij sin N j
K (ϕi − ϕj ) 2
− funknown (ϕi )KS sin (Kϕi ) X ′ + funknown (ϕi )KG Wij g(ϕi , ϕj ),
(13)
j ′
with KG being a scaling factor, Wij a weight matrix that does not differentiate between knowns and unknowns, i.e., ( ′
Wij =
−1, i ̸= j, i ∧ j ∈ row, column or box, 0, else,
(14)
and g(ϕi , ϕj ) producing a non-vanishing value if and only if ϕi and ϕj are the n multiples of 2π far apart from each (i.e. they output the same digit) or in mathematical terms: ( 1, |ϕi − ϕj | ≈ 2πn, g(ϕi , ϕj ) = (15) 0, else. We emulate g(ϕi , ϕj ) with the following function: (1 − cos (ϕi − ϕj )) , g(ϕi , ϕj ) = exp − σ
(16)
where σ tunes the width of the peaks around the n-th multiples of 2π. B. Benchmark set-up The model presented in Eq. 13 has four tuneable parameters, namely: • the weight between two unknowns wu , • the scaling factor of the K-th harmonics KS , • the scaling factor of the g(ϕi , ϕj ) function KG , • the width of the g(ϕi , ϕj ) function’s peaks σ. Given the amount of tuneable parameters, a thorough sweep must first be conducted to identify the appropriate combinations of parameters. Hence, we generate a training set with a size of 250 samples to narrow down the space of parameter combinations and finally utilize a test set with a size of 1000 samples to report on the best achieved performance. The Sudoku puzzles are created with the PySudoku library [39]. In particular, we benchmark our proposed ONN-Suduku solver (see Eq. 13) on Sudoku puzzles of sizes 4 × 4 and 9 × 9 for seven different unknown ratios (i.e. how many of the digits in the puzzle are not known) [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] against an HNN Sudoku-solver [36] and an already established ONN
Sudoku-solver [28] in terms of accuracy. Finally, we report the order parameters recorded for each unknown ratio and each size to assess how well the model converged towards the preestablished roots. All benchmarks are performed in Python. IV. R ESULTS First, we start with the accuracy benchmarks for 4 × 4 and 9 × 9 Sudoku puzzles. In particular, we report on the best accuracies achieved with our ONN- and the HNN-Sudoku solver [36] for seven different unknown ratios [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] on 1000 Sudoku puzzles. Additionally, we include the accuracies of the already developed ONN-Sudoku solver for five different unknown ratios [0.1, 0.2, 0.3, 0.4, 0.5] from the authors in [28]. Figures 3a and 3b summarize all these accuracies as a function of the unknown ratio for sizes 4 × 4 and 9 × 9, respectively. Let us start by analyzing the results for the size 4 × 4. Regarding the HNN-Sudoku solver, as the unknown ratio increases, performance rapidly deteriorates, effectively being unable to solve any puzzles correctly for unknown ratios greater than 60%. The accuracy drastically improves with the conventional ONN-Sudoku solver achieving almost 80% for an unknown ratio of 50%. In contrast, our proposed approach significantly outperforms both established methods, since it remains almost flawless throughout the benchmarks with the accuracy slightly decreasing for an unknown ratio of 90%. We speculate that this vastly superior performance can be traced back to the additional term presented in Eq. 13 as it ensures a steady state if and only if all cells are correctly filled out. Let us now move on to the benchmarks on 9×9 Sudoku puzzles. Fully analogously to previously reported benchmarks, the accuracy of the HNN solver quickly worsens as the number of unknown digits grows, effectively being incapable of solving any puzzles correctly for unknown ratios greater than 40%. While the conventional ONN-Sudoku solver starts off with an almost perfect accuracy, the accuracy rapidly decreases, reporting no solved puzzles for an unknown ratio of 50%. As far as our proposed ONN-Sudoku solver is concerned, while still producing superior accuracies in comparison to the two other methods, it is no longer able to maintain a flawless accuracy as the number of unknowns ramps up. In particular, although perfect accuracies are reached for unknown ratios up to 25% and an accuracy of more than 80% is obtained for an unknown ratio of 37.5%, the precision starts to decrease, achieving approximately 50% for an unknown ratio of 50% and reaching the limits of the model at unknown ratios greater than 70% by not being able to solve any puzzles correctly. We assume that the order parameters will decrease as the number of unknowns increases, implying that either the additional term presented in Eq. 13 has not a strong enough influence and/or the model needs more oscillatory cycles to reach a stable state. In the second part of the benchmarks, let us closely examine the order parameters of our proposed ONN-Sudoku solver for 4 × 4 and 9 × 9 Sudoku puzzles. In particular, we report on the order parameters associated with the best accuracies reached by our solver for seven different unknown
1.0
0.8
0.8
0.6
ONN - proposed approach ONN - conventional approach HNN
0.4 0.2 0.0
Accuracy
Accuracy
1.0
ONN - proposed approach ONN - conventional approach HNN
0.6 0.4 0.2
0.2
0.4
0.6
Unknown ratio
0.8
(a)
0.0
0.2
0.4
0.6
Unknown ratio
0.8
(b)
Fig. 3: Accuracy as a function of unknown ratio of the HNN-Sudoku solver [36], the conventional ONN-Sudoku solver [28] and our proposed ONN-Sudoku solver for (a) 4×4 and (b) 9×9 Sudoku puzzles, respectively. For our ONN- and the HNN-Sudoku solver, we present the accuracies achieved on 1000 samples and for the conventional ONN approach, we plot the accuracy values reported in [28].
ratios [0.125, 0.25, 0.375, 0.5, 0.625, 0.7, 0.875] on 1000 Sudoku puzzles. Figures 4a and 4b summarize all these order parameters as histograms for sizes 4×4 and 9×9, respectively. Clearly, for 4 × 4 Sudokus, where almost perfect accuracies are achieved, the oscillators are able to converge consistently to their pre-established roots for all seven different unknown ratios. In contrast, for 9 × 9 Sudokus, as the number of unknown digits ramps up, the distribution of order parameters starts shifting towards lower values. In particular, while for an unknown ratio of 0.125, almost all order parameters are clustered around 1, for an unknown ratio of 0.875, the vast majority of order parameters is below 0.5. As theorized earlier, a high accuracy induces a high order parameter, while lower accuracies are associated with lower order parameters. This relationship can be attributed to the g(ϕi , ϕj ) function in Eq. 13, as it ensures a stable state if and only if no Sudoku rules are violated. A low order parameter may hint at either the additional term not being scaled appropriately and/or the ONN requiring more oscillatory cycles to reach a stable state. V. D ISCUSSION In this paper, we present our developed ONN-Sudoku solver inspired by the Graph Coloring problem. In Python, we successfully demonstrated that an appropriate combination of tuneable parameters can lead to almost flawless precisions for 4 × 4 Sudokus and significantly improved accuracies for 9 × 9 Sudokus. However, it needs to be acknowledged that the model we considered (see Eq. 13) has a few drawbacks. In particular, the model presented in this paper has four tuneable parameters, thus narrowing the search down through a parameter sweep to a suitable combination of parameters can be quite costly and time intensive. Due to hardware restrictions and the parameter sweep, we opted for a training set size of only 250 samples.
Hence, future versions of this Sudoku solver should look into eliminating a few of the tuneable parameters. Furthermore, while our proposed Sudoku solver outperformed the existing HNN- and ONN-solvers, it started hitting its limits for 9 × 9 and unknown ratios larger than 50%. As the order parameter histograms confirmed, there is a clear relationship between accuracy and order parameter, as high accuracies induce high order parameters and low accuracies occur only with low order parameters. Based on these results, we speculate that the additional term with the g(ϕi , ϕj ) function (see Eq. 13) is not scaled properly and/or the ONN did not have enough oscillatory cycles to settle. Additionally, it may be worth looking into alternative g(ϕi , ϕj ) function formulations than the one presented in Eq. 16. While our model deploys simple, sinusoidal KuramotoONNs [14, 37], there has recently been a push in the community to harness the dynamics of highly nonlinear oscillators for computing [40–42]. Although it is not a limitation per se, it would be worth investigating, how a highly nonlinear or even noisy oscillator would impact the accuracy as well as the order parameter compared to the noiseless Kuramoto model deployed here [13, 43–47]. VI. C ONCLUSIONS In this paper, we propose a Sudoku solver based on Oscillatory Neural Networks (ONNs) by formulating the problem as a Graph Coloring problem. Sudokus embody a well-established logical puzzle, which can be seen as a constrained combinatorial optimization problem (COP). In particular, as the Sudoku rules forbid a placement of identical digits in common rows, columns, and boxes, Graph Coloring, a problem which boils down to assigning “colors” to graph nodes such that no two
0.6
800 600 400 200 0
0.8
0
1.0
0.2
0.4
0.6
1000 750 500 250 0
0.4
0.2
0.4
0.6
0.8
1.0
Order parameter
Count
150
0.8
1.0
1000 750 500 250 0
0.2
0.4
Unknown ratio = 0.625
50 0.0
0.2
0.4
0.6
0.8
Order parameter
1.0
0.0
0.2
0.0
0.6
0.2
0.4
0.6
0.8
0.6
800 600 400 200 0
(a)
1.0
800 600 400 200 0
0.0
0.2
0.6
0.0
0.2
0.8
Order parameter
0.4
0.6
0.8
Order parameter
1.0
0.0
0.2
0.4
0.6
0.8
0.4
0.6
100 75 50 25 0
1.0
Order parameter
Unknown ratio = 0.5
0.8
400 200 0
1.0
0.0
0.2
0.4
0.6
0.8
Order parameter
1.0
Unknown ratio = 0.875
50 0.4
250 0
1.0
Order parameter
100
0.2
500
Unknown ratio = 0.375
Unknown ratio = 0.75
0.0
0.8
750
Unknown ratio = 0.875
1.0
Order parameter
0.8
150
0
0.4
Order parameter
Unknown ratio = 0.75
Order parameter
100
0
1.0
Unknown ratio = 0.25
0.0
200 0
Count
Unknown ratio = 0.125
0.0
0.6
Order parameter
Count
1000 750 500 250 0
0.2
Count
Count
Unknown ratio = 0.625
0.0
0.8
Order parameter
Count
Count 0.0
400
Unknown ratio = 0.5
1000
Count
0.4
Order parameter
200
Unknown ratio = 0.375
Count
0.2
400
Count
0.0
600
Count
Unknown ratio = 0.25 Count
Unknown ratio = 0.125
Count
Count
1000 750 500 250 0
1.0
0.0
0.2
0.4
0.6
0.8
Order parameter
1.0
(b)
Fig. 4: Order parameters associated with the best accuracies achieved by our proposed ONN-Sudoku solver plotted as histograms for seven different unknown ratios [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875] as well as (a) 4 × 4 and (b) 9 × 9 Sudoku puzzles, respectively.
adjacent nodes possess the same color and the number of unique colors is minimized, presents a natural COP fit. Oscillatory Neural Networks (ONNs) present an attractive physics-based paradigm for solving hard combinatorial optimization problems by harnessing the rich dynamics of coupled oscillators. While there exists a quite computationally expensive Graph Coloring mapping to ONNs, it cannot simply be utilized for Sudokus as it allows potential rule violations. Instead, our ONN-Sudoku solver expands a modified, computationally cheaper version of Graph Coloring by an additional term which ensures that no two oscillators sharing a row, column or box yield identical digits. We benchmark our ONN-Sudoku solver for seven different numbers of unknown digits on 4 × 4 and 9 × 9 puzzles in terms of accuracy and order parameter against two solvers a solver based on Hopfield Neural Networks and an already established ONN-solver. The accuracy measurements prove that our proposed solver significantly outperforms both existent solvers on 4 × 4 as well as 9 × 9 puzzles, yielding almost flawless accuracies for 4 × 4 Sudokus, while hitting a limit on 9 × 9 Sudokus for large unknown ratios. The order
parameter benchmarks reveal that high accuracies are associated with large order parameter values, whereas degrading precisions occur with lower order parameters, hinting at either the additional term not being scaled properly and/or the model requiring more oscillatory cycles to stabilize. This paper demonstrates the relevance of expanding the ODE by an additional term in ONNs for constrained combinatorial optimization to kick the oscillators out of undesired states and produce vanishing gradients for only configurations satisfying all the constraints. In particular, this term allowed ONNs to outperform already established Sudoku solvers based on ONNs and HNNs in terms of accuracy. Last but certainly not least, to improve the performance, future work should investigate the benefits of highly nonlinear oscillators and different additional term function formulations. R EFERENCES [1] X.-S. Yang and S. Koziel, Computational optimization and applications in engineering and industry. Springer Science & Business Media, 2011, vol. 359.
[2] A. R. Yıldız, “An effective hybrid immune-hill climbing optimization approach for solving design and manufacturing optimization problems in industry,” Journal of Materials Processing Technology, vol. 209, no. 6, pp. 2773–2780, 2009. [3] M. Fathi, M. Khakifirooz, and P. M. Pardalos, Optimization in large scale problems: Industry 4.0 and Society 5.0 applications. Springer, 2019, vol. 152. [4] P. Bangert, Optimization for industrial problems. Springer Science & Business Media, 2012. [5] R. M. Karp, “Reducibility among combinatorial problems,” in 50 Years of Integer Programming 1958-2008: from the Early Years to the State-of-the-Art. Springer, 2009, pp. 219–241. [6] J. von Neumann, “First draft of a report on the edvac,” IEEE Annals of the History of Computing, vol. 15, no. 4, pp. 27–75, 1993. [7] D. Efnusheva, A. Cholakoska, and A. Tentov, “A survey of different approaches for overcoming the processormemory bottleneck,” International Journal of Computer Science and Information Technology, vol. 9, no. 2, pp. 151–163, 2017. [8] R. Toews, “Deep learning’s carbon emissions problem,” Forbes. [Online]. Available: https://www.forbes.com/sites/robtoews/2020/06/17/deeplearnings-climate-change-problem/ [9] A. Todri-Sanial, C. Delacour, M. Abernot, and F. Sabo, “Computing with oscillators from theoretical underpinnings to applications and demonstrators,” Npj unconventional computing, vol. 1, no. 1, p. 14, 2024. [10] F. Rosenblatt, “The perceptron: a probabilistic model for information storage and organization in the brain.” Psychological review, vol. 65, no. 6, p. 386, 1958. [11] B. Yegnanarayana, Artificial neural networks. PHI Learning Pvt. Ltd., 2009. [12] J. Zou, Y. Han, and S.-S. So, “Overview of artificial neural networks,” Artificial neural networks: methods and applications, pp. 14–22, 2009. [13] E. M. Izhikevich, “Simple model of spiking neurons,” IEEE Transactions on neural networks, vol. 14, no. 6, pp. 1569–1572, 2003. [14] T. Wang and J. Roychowdhury, “Oim: Oscillator-based ising machines for solving combinatorial optimisation problems,” in Unconventional Computation and Natural Computation: 18th International Conference, UCNC 2019, Tokyo, Japan, June 3–7, 2019, Proceedings 18. Springer, 2019, pp. 232–256. [15] J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proceedings of the National Academy of Sciences of the United States of America, vol. 79, pp. 2554–8, 05 1982. [16] H. Ramsauer et al., “Hopfield networks is all you need,” 2021. [17] E. Ising, “Beitrag zur theorie des ferro-und paramagnetismus,” Ph.D. dissertation, Grefe & Tiedemann Hamburg, Germany, 1924.
[18] A. Lucas, “Ising formulations of many np problems,” Frontiers in physics, vol. 2, p. 74887, 2014. [19] F. C. Hoppensteadt and E. M. Izhikevich, “Pattern recognition via synchronization in phase-locked loop neural networks,” IEEE Transactions on Neural Networks, vol. 11, no. 3, pp. 734–738, 2000. [20] B. F. Haverkort and A. Todri-Sanial, “Overcoming quadratic hardware scaling for a fully connected digital oscillatory neural network,” Frontiers in Neuroscience, vol. Volume 19 - 2025, 2026. [Online]. Available: https://www.frontiersin.org/journals/neuroscience/articles/10.3389/fnin [21] M. Abernot, T. Gil, M. Jiménez, J. Núñez, M. J. Avellido, B. Linares-Barranco, T. Gonos, T. Hardelin, and A. Todri-Sanial, “Digital implementation of oscillatory neural network for image recognition applications,” Frontiers in Neuroscience, vol. 15, p. 713054, 2021. [22] O. Maher, F. Horst, F. Sabo, D. Fancheng, V. Bragaglia, S. Karg, M. Sousa, A. Todri-Sanial, and B. J. Offrein, “Implementation of fitzhugh-nagumo neurons using nanoscale vo 2 devices,” IEEE, pp. 729–732, 2024. [23] A. Gower, “Learning at the speed of physics: Equilibrium propagation on oscillator ising machines,” 2025. [Online]. Available: https://arxiv.org/abs/2510.12934 [24] F. Sabo and A. Todri-Sanial, “Classonn: Classification with oscillatory neural networks using the kuramoto model,” IEEE, pp. 1–2, 2024. [25] M. Abernot and A. Todri-Sanial, “Training energy-based single-layer hopfield and oscillatory networks with unsupervised and supervised algorithms for image classification,” Neural Computing and Applications, vol. 35, no. 25, pp. 18 505–18 518, 2023. [26] N. R. Rohan, C. Vigneswaran, S. Ghosh, K. Rajendran, A. Gaurav, and V. S. Chakravarthy, “Deep oscillatory neural network,” Scientific Reports, vol. 15, no. 1, p. 40968, 2025. [Online]. Available: https://doi.org/10.1038/s41598-025-24837-4 [27] T. Miyato, S. Löwe, A. Geiger, and M. Welling, “Artificial kuramoto oscillatory neurons,” 2025. [Online]. Available: https://arxiv.org/abs/2410.13821 [28] B. F. Haverkort, F. Sbravati, S. Porfir, and A. TodriSanial, “Solving sudoku using oscillatory neural networks,” Neuromorphic Computing and Engineering, 2025. [29] L. Sun, M. X. Burns, and M. C. Huang, “General oscillator-based Ising-machine models with phase-amplitude dynamics and polynomial interactions,” Physical Review Applied, vol. 24, no. 4, p. 044035, Oct. 2025. [Online]. Available: https://link.aps.org/doi/10.1103/7s4v-4hs4 [30] Y. E. Gonul, C. E. Kayan, I. Mustafazade, N. Kandasamy, and B. Taskin, “Gpu-accelerated simulated oscillator ising/potts machine solving combinatorial optimization problems,” 2025. [Online]. Available: https://arxiv.org/abs/2505.22631 [31] S. K. Vadlamani, T. P. Xiao, and E. Yablonovitch, “Combinatorial optimization using the Lagrange
primal-dual dynamics of parametric oscillator & Sons, 2011. networks,” Physical Review Applied, vol. 21, [47] E. Panteley, A. Loria, and A. El Ati, “On the stabilno. 4, p. 044042, Apr. 2024. [Online]. Available: ity and robustness of stuart-landau oscillators,” IFAChttps://link.aps.org/doi/10.1103/PhysRevApplied.21.044042 PapersOnLine, vol. 48, no. 11, pp. 645–650, 2015. [32] M. K. Bashar, Z. Li, V. Narayanan, and N. Shukla, “An fpga-based max-k-cut accelerator exploiting oscillator synchronization model,” IEEE, pp. 1–8, 2024. [33] O. Maher, M. Jiménez, C. Delacour, N. Harnack, J. Núñez, M. J. Avedillo, B. Linares-Barranco, A. Todri-Sanial, G. Indiveri, and S. Karg, “A CMOS-compatible oscillation-based VO2 Ising machine solver,” Nature Communications, vol. 15, no. 1, p. 3334, Apr. 2024. [Online]. Available: https://www.nature.com/articles/s41467-024-47642-5 [34] T. Kanao and H. Goto, “Simulated bifurcation assisted by thermal fluctuation,” Communications Physics, vol. 5, no. 1, p. 153, Jun. 2022. [Online]. Available: https://www.nature.com/articles/s42005-022-00929-9 [35] A. Mallick, M. K. Bashar, Z. Lin, and N. Shukla, “Computational models based on synchronized oscillators for solving combinatorial optimization problems,” Physical Review Applied, vol. 17, no. 6, p. 064064, 2022. [36] V. Mladenov, P. Karampelas, C. Pavlatos, and E. Zirintsis, “Solving sudoku puzzles by using hopfield neural networks,” Proc. of ICACM, vol. 11, pp. 174–179, 2011. [37] Y. Kuramoto, Chemical turbulence. Springer, 1984. [38] Y. Cheng, M. Khairul Bashar, N. Shukla, and Z. Lin, “A control theoretic analysis of oscillator ising machines,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 34, no. 7, 2024. [39] J. Sieu, “py-sudoku,” https://github.com/jeffsieu/pysudoku, 2024. [40] F. Sabo, F. Du, N. Dinç, and A. Todri-Sanial, “Harnessing dynamics of van der pol oscillators for phase computing with oscillatory neural networks,” 2026. [41] M. R. E. U. Shougat, X. Li, T. Mollik, and E. Perkins, “An information theoretic study of a duffing oscillator array reservoir computer,” Journal of Computational and Nonlinear Dynamics, vol. 16, no. 8, p. 081004, 2021. [42] O. V. Maslennikov, D. S. Shchapin, and V. I. Nekorkin, “Binary classification via spatiotemporal dynamics in reservoir computing rings of fitzhugh–nagumo neurons,” The European Physical Journal Special Topics, vol. 234, no. 15, pp. 4033–4041, 2025. [43] B. van der Pol, “Lxxxviii. on “relaxation-oscillations”,” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, vol. 2, no. 11, pp. 978–992, 1926. [44] R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical journal, vol. 1, no. 6, pp. 445–466, 1961. [45] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proceedings of the IRE, vol. 50, no. 10, pp. 2061–2070, 1962. [46] I. Kovacic and M. J. Brennan, The Duffing equation: nonlinear oscillators and their behaviour. John Wiley