Graph Coloring Approach to Solving Sudoku with Oscillatory Neural Networks

# Graph Coloring Approach to Solving Sudoku with Oscillatory Neural Networks ††thanks: This 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.
Source: [https://arxiv.org/html/2607.15814](https://arxiv.org/html/2607.15814)
###### 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 HNN\- and ONN solvers in terms of accuracy\. In particular, we are able to achieve nearly flawless accuracies on4×44\\times 4as well as rather high accuracies on9×99\\times 9Sudoku puzzles for different numbers of unknown digits\.

## IIntroduction

Constrained combinatorial optimization problems \(COPs\) are ubiquitous in the industry, ranging from scheduling through chip design to the finance sector\[[44](https://arxiv.org/html/2607.15814#bib.bib30),[46](https://arxiv.org/html/2607.15814#bib.bib27),[7](https://arxiv.org/html/2607.15814#bib.bib28),[3](https://arxiv.org/html/2607.15814#bib.bib29)\]\. 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\[[18](https://arxiv.org/html/2607.15814#bib.bib1)\]\. 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\[[42](https://arxiv.org/html/2607.15814#bib.bib35),[6](https://arxiv.org/html/2607.15814#bib.bib36),[39](https://arxiv.org/html/2607.15814#bib.bib37)\]\. 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\)\[[38](https://arxiv.org/html/2607.15814#bib.bib7)\]\.

In contrast to conventional Artificial Neural Networks\[[32](https://arxiv.org/html/2607.15814#bib.bib38),[45](https://arxiv.org/html/2607.15814#bib.bib40),[47](https://arxiv.org/html/2607.15814#bib.bib39)\], Oscillatory Neural Networks present an emerging physics\-based computing paradigm rooted in the dynamics of a network of typically fully coupled oscillators aiming to minimize an underlying energy function\[[16](https://arxiv.org/html/2607.15814#bib.bib18),[43](https://arxiv.org/html/2607.15814#bib.bib5),[38](https://arxiv.org/html/2607.15814#bib.bib7)\]\. In such networks, the information is encoded in the phases of the oscillators\. Thanks to their close ties to Hopfield Neural Networks\[[13](https://arxiv.org/html/2607.15814#bib.bib42),[30](https://arxiv.org/html/2607.15814#bib.bib43)\]and the Ising model\[[15](https://arxiv.org/html/2607.15814#bib.bib41),[21](https://arxiv.org/html/2607.15814#bib.bib2)\], ONNs are well suited for auto\-associative memory tasks\[[14](https://arxiv.org/html/2607.15814#bib.bib4),[11](https://arxiv.org/html/2607.15814#bib.bib44),[1](https://arxiv.org/html/2607.15814#bib.bib45),[22](https://arxiv.org/html/2607.15814#bib.bib46),[10](https://arxiv.org/html/2607.15814#bib.bib47),[34](https://arxiv.org/html/2607.15814#bib.bib48),[2](https://arxiv.org/html/2607.15814#bib.bib49),[31](https://arxiv.org/html/2607.15814#bib.bib50),[26](https://arxiv.org/html/2607.15814#bib.bib51)\]and combinatorial optimization problems\[[43](https://arxiv.org/html/2607.15814#bib.bib5),[12](https://arxiv.org/html/2607.15814#bib.bib11),[37](https://arxiv.org/html/2607.15814#bib.bib52),[9](https://arxiv.org/html/2607.15814#bib.bib53),[40](https://arxiv.org/html/2607.15814#bib.bib54),[4](https://arxiv.org/html/2607.15814#bib.bib22),[23](https://arxiv.org/html/2607.15814#bib.bib25),[17](https://arxiv.org/html/2607.15814#bib.bib55),[24](https://arxiv.org/html/2607.15814#bib.bib3)\]\. In this paper, we propose an ONN\-based 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 sizeN×NN\\times Nand some known digits, the aim lies in filling out the remaining empty cells of the grid with integers from the range\[1,N\]\[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\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]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 for4×44\\times 4,9×99\\times 9as well as16×1616\\times 16puzzles in terms of accuracy\[[27](https://arxiv.org/html/2607.15814#bib.bib10)\]\. In contrast, in this paper, we present an alternative ONN\-Sudoku 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 on4×44\\times 4and9×99\\times 9Sudoku puzzles by reaching almost perfect accuracies on4×44\\times 4and predominantly high accuracies on9×99\\times 9Sudokus for various numbers of unknown digits\.

This paper is structured as follows\. In Section[II](https://arxiv.org/html/2607.15814#S2), we introduce the key concepts behind physics\-based computing with Oscillatory Neural Networks \(ONNs\)\. Furthermore, Section[III](https://arxiv.org/html/2607.15814#S3)outlines our proposed ONN\-Sudoku solver based on Graph Coloring together with the benchmarking set\-up\. Most importantly, our benchmarks on4×44\\times 4and9×99\\times 9Sudokus in terms of accuracy against the established HNN\-\[[27](https://arxiv.org/html/2607.15814#bib.bib10)\]and ONN\-Sudoku solvers\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]are presented in Section[IV](https://arxiv.org/html/2607.15814#S4)\. Last but certainly not least, the results are discussed in Section[V](https://arxiv.org/html/2607.15814#S5)and the paper is concluded with Section[VI](https://arxiv.org/html/2607.15814#S6)\.

## IIBackground

### II\-AIsing model

This paper revolves around Oscillatory Neural Networks \(ONNs\), an emerging physics\-based paradigm rooted in the dynamics of coupled oscillators\[[38](https://arxiv.org/html/2607.15814#bib.bib7)\]\. In contrast to Artificial Neural Networks\[[32](https://arxiv.org/html/2607.15814#bib.bib38),[45](https://arxiv.org/html/2607.15814#bib.bib40),[47](https://arxiv.org/html/2607.15814#bib.bib39)\], ONN computing boils down to finding global minima in an underlying energy function\[[16](https://arxiv.org/html/2607.15814#bib.bib18),[43](https://arxiv.org/html/2607.15814#bib.bib5),[38](https://arxiv.org/html/2607.15814#bib.bib7)\]\. The ONNs are inherently linked to the Ising model\[[15](https://arxiv.org/html/2607.15814#bib.bib41),[21](https://arxiv.org/html/2607.15814#bib.bib2)\], a mathematical model of ferromagnetism describing a spin configuration of a system\[[15](https://arxiv.org/html/2607.15814#bib.bib41)\]\. In particular, the dynamics of the model can be described with the following Hamiltonian:

H=−∑\(i,j\)Ji​j​σi​σj\.H=\-\\sum\_\{\(i,j\)\}J\_\{ij\}\\sigma\_\{i\}\\sigma\_\{j\}\.\(1\)Given a coupling weight matrixJi​jJ\_\{ij\}, the goal lies in identifying the natural state of the system by minimizing the outlined Hamiltonian through the set of binary spinsσi\\sigma\_\{i\}\.

### II\-BOscillatory Neural Networks \(ONNs\)

One can bridge the gap between the Ising model and ONNs by mapping the spinsσi∈\{1,−1\}\\sigma\_\{i\}\\in\\\{1,\-1\\\}to the oscillators’ phasesϕi∈\{0,π\}\\phi\_\{i\}\\in\\\{0,\\pi\\\}resulting in the corresponding ONN\-Hamiltonian\. Moreover, the ordinary differential equation \(ODE\) describing the dynamics ofNNoscillators in an ONN can thus be determined simply by following the negative gradient of the ONN\-Hamiltonian\[[14](https://arxiv.org/html/2607.15814#bib.bib4),[43](https://arxiv.org/html/2607.15814#bib.bib5)\]:

ϕi˙=1N​∑jWi​j​sin⁡\(\(ϕj−ϕi\)\),\\dot\{\\phi\_\{i\}\}=\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\(\\phi\_\{j\}\-\\phi\_\{i\}\)\},\(2\)whereWi​jW\_\{ij\}stands for the coupling between theii\-th andjj\-th oscillator\. The ODE in Eq\.[2](https://arxiv.org/html/2607.15814#S2.E2)is often called the Kuramoto model\[[20](https://arxiv.org/html/2607.15814#bib.bib6)\]\. As already mentioned, the phases \(0,π0,\\pi\) encode the Booleans \(1,−11,\-1\)\. However, when running an ONN, it rarely happens that the phases settle at the pre\-established binary values0andπ\\pi, but rather at continuous ones\. Hence, to force the oscillators into the binary phases, one introduces an additional term scaled by a factorKSK\_\{S\}into the ODE called*Subharmonic Injection Locking*\(SHIL\) which injects the second harmonics into the oscillators\[[38](https://arxiv.org/html/2607.15814#bib.bib7),[5](https://arxiv.org/html/2607.15814#bib.bib8),[43](https://arxiv.org/html/2607.15814#bib.bib5)\]:

ϕi˙=1N​∑jWi​j​sin⁡\(\(ϕj−ϕi\)\)−KS​sin⁡\(\(2​ϕi\)\)\.\\dot\{\\phi\_\{i\}\}=\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\(\\phi\_\{j\}\-\\phi\_\{i\}\)\}\-K\_\{S\}\\sin\{\(2\\phi\_\{i\}\)\}\.\(3\)

## IIIMethods

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\.

### III\-AProposed ONN\-Sudoku solver

#### III\-A1Graph 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\[[18](https://arxiv.org/html/2607.15814#bib.bib1)\], Graph Coloring is a well\-known hard combinatorial optimization problem\. As shown by the author in\[[21](https://arxiv.org/html/2607.15814#bib.bib2)\], all NP\-hard problems can be mapped to a corresponding Ising Hamiltonian with Graph Coloring not being an exception\.

Furthermore, authors in\[[24](https://arxiv.org/html/2607.15814#bib.bib3)\]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\-KK\-Cut with the minimalKKsuch that no two adjacent nodes share the same set\. In the case of a Sudoku, the number of sets orKKis well\-known and does not need to be minimized, as it corresponds to the size of the Sudoku\. Thus, the desired phasesϕideal\\phi\_\{\\text\{ideal\}\}for a sizeKKSudoku read as follows:

ϕideal=\{2​π​\(k−1\)K\|k∈\{1,2,…,K\}\}\.\\phi\_\{\\text\{ideal\}\}=\\biggl\\\{\\frac\{2\\pi\(k\-1\)\}\{K\}\\bigg\|k\\in\\\{1,2,\.\.\.,K\\\}\\biggr\\\}\.\(4\)
Hence, for a9×99\\times 9Sudoku \(K=9K=9\), we get the following mapping of the numbers to the phases:

1→02→2​π93→4​π94→2​π35→8​π96→10​π97→4​π38→14​π99→16​π9\\begin\{split\}1&\\rightarrow 0\\\\ 2&\\rightarrow\\frac\{2\\pi\}\{9\}\\\\ 3&\\rightarrow\\frac\{4\\pi\}\{9\}\\\\ 4&\\rightarrow\\frac\{2\\pi\}\{3\}\\\\ 5&\\rightarrow\\frac\{8\\pi\}\{9\}\\\\ 6&\\rightarrow\\frac\{10\\pi\}\{9\}\\\\ 7&\\rightarrow\\frac\{4\\pi\}\{3\}\\\\ 8&\\rightarrow\\frac\{14\\pi\}\{9\}\\\\ 9&\\rightarrow\\frac\{16\\pi\}\{9\}\\end\{split\}\(5\)
Let us now turn our attention towards Max\-KK\-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\-KK\-Cut, a partition intoKKsets with the maximum sum of the cut edges shall be found\. The authors in\[[24](https://arxiv.org/html/2607.15814#bib.bib3)\]solve Max\-KK\-Cut with Oscillatory Neural Networks by extending ODE from Eq\.[3](https://arxiv.org/html/2607.15814#S2.E3)to the following form:

ϕi˙=−1N​∑jWi​j​sin⁡\(\(ϕi−ϕj\+f​\(ϕi−ϕj\)\)\)−KS​sin⁡\(\(K​ϕi\)\),\\begin\{split\}\\dot\{\\phi\_\{i\}\}=&\-\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\(\\phi\_\{i\}\-\\phi\_\{j\}\+f\(\\phi\_\{i\}\-\\phi\_\{j\}\)\)\}\\\\ &\-K\_\{S\}\\sin\{\(K\\phi\_\{i\}\)\},\\end\{split\}\(6\)withf​\(Δ​ϕi​j\)f\(\\Delta\\phi\_\{ij\}\)ensuring that the ideal phases from Eq\.[4](https://arxiv.org/html/2607.15814#S3.E4)produce a vanishing gradient when acted upon with a sine function\. Furthermore, the functionffreads as follows:

f​\(Δ​ϕi​j\)=∑k=1K−1\(\(2​k−1\)​π−2​k​πK\)\(exp⁡\(−\(Δ​ϕi​j−2​k​πK\)22​σ2\)−exp⁡\(−\(Δ​ϕi​j\+2​k​πK\)22​σ2\)\),\\begin\{split\}f\(\\Delta\\phi\_\{ij\}\)=&\\sum\_\{k=1\}^\{K\-1\}\\left\(\(2k\-1\)\\pi\-\\frac\{2k\\pi\}\{K\}\\right\)\\\\ &\\biggl\(\\exp\{\-\\frac\{\(\\Delta\\phi\_\{ij\}\-\\frac\{2k\\pi\}\{K\}\)^\{2\}\}\{2\\sigma^\{2\}\}\}\\\\ &\-\\exp\{\-\\frac\{\(\\Delta\\phi\_\{ij\}\+\\frac\{2k\\pi\}\{K\}\)^\{2\}\}\{2\\sigma^\{2\}\}\}\\biggl\),\\end\{split\}\(7\)withσ\\sigmabeing a tuneable parameter\. Fig\.[1](https://arxiv.org/html/2607.15814#S3.F1)depicts the functionf​\(Δ​ϕi​j\)f\(\\Delta\\phi\_\{ij\}\)from Eq\.[7](https://arxiv.org/html/2607.15814#S3.E7)and the termsin⁡\(\(Δ​ϕi​j\+f​\(Δ​ϕi​j\)\)\)\\sin\{\(\\Delta\\phi\_\{ij\}\+f\(\\Delta\\phi\_\{ij\}\)\)\}forK=4K=4and a relatively largeσ\\sigma\(=0\.15=0\.15\) to showcase the form of the function\.

![Refer to caption](https://arxiv.org/html/2607.15814v1/x1.png)Figure 1:Functionf​\(Δ​ϕi​j\)f\(\\Delta\\phi\_\{ij\}\)from Eq\.[7](https://arxiv.org/html/2607.15814#S3.E7)\(left\) and termsin⁡\(\(Δ​ϕi​j\+f​\(Δ​ϕi​j\)\)\)\\sin\{\(\\Delta\\phi\_\{ij\}\+f\(\\Delta\\phi\_\{ij\}\)\)\}\(right\) forK=4K=4\.Clearly, utilizing the model from Eq\.[6](https://arxiv.org/html/2607.15814#S3.E6)has multiple drawbacks\. First of all, the function in Eq\.[7](https://arxiv.org/html/2607.15814#S3.E7)has a tuneable parameterσ\\sigma, which controls the width of the peaks\. Second, computing the terms in Eq\.[7](https://arxiv.org/html/2607.15814#S3.E7)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 termsin⁡\(\(Δ​ϕi​j\+f​\(Δ​ϕi​j\)\)\)\\sin\{\(\\Delta\\phi\_\{ij\}\+f\(\\Delta\\phi\_\{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:

ϕi˙=−1N​∑jWi​j​sin⁡\(\(K2​\(ϕi−ϕj\)\)\)−KS​sin⁡\(\(K​ϕi\)\)\.\\dot\{\\phi\_\{i\}\}=\-\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\\left\(\\frac\{K\}\{2\}\(\\phi\_\{i\}\-\\phi\_\{j\}\)\\right\)\}\-K\_\{S\}\\sin\{\(K\\phi\_\{i\}\)\}\.\(8\)In particular, the first term ensures that the oscillators are2​π​nK\\frac\{2\\pi n\}\{K\}apart, while the second term guarantees that the oscillators settle at2​π​nK\\frac\{2\\pi n\}\{K\}values\. We utilize Eq\.[8](https://arxiv.org/html/2607.15814#S3.E8)instead of Eq\.[6](https://arxiv.org/html/2607.15814#S3.E6)in our Sudoku approach\.

Based on Eq\.[8](https://arxiv.org/html/2607.15814#S3.E8), we define an order parameterκ\\kappato track how well the oscillators converged to the desired phases as outlined in Eq\.[4](https://arxiv.org/html/2607.15814#S3.E4):

κ=1N​∑i1−\|Im​\(exp⁡\(\(i​K​ϕi2\)\)\)\|\.\\kappa=\\frac\{1\}\{N\}\\sum\_\{i\}1\-\\left\|\\text\{Im\}\\left\(\\exp\{\\left\(\\frac\{iK\\phi\_\{i\}\}\{2\}\\right\)\}\\right\)\\right\|\.\(9\)
The order parameter is bounded to the range\[0,1\]\[0,1\], where values close to11indicate that the phasesϕi\\phi\_\{i\}converged well to the pre\-established phasesϕi​d​e​a​l\\phi\_\{ideal\}, whereas values near0hint at a poor convergence toϕi​d​e​a​l\\phi\_\{ideal\}\. Let us now discuss the form of the weight matrixWi​jW\_\{ij\}\.

#### III\-A2Weight 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\-KK\-Cut with the minimalKKsuch 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\-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 itwuw\_\{u\}\) shall be negative and have a smaller magnitude than the weight connecting a known \(−1\-1in our case\)\. Putting all of these conditions together, we arrive at the following form of the weight matrixWi​jW\_\{ij\}:

Wi​j=\{−1,i≠j, i∧j∈row, column or box, i∨j∈known,wu,i≠j, i∧j∈row, column or box, i∧j∉known,0,else,W\_\{ij\}=\\begin\{cases\}\-1,\\text\{i $\\neq$ j, i $\\land$ j $\\in$ row, column or box, i $\\lor$ j $\\in$ known,\}\\\\ w\_\{u\},\\text\{i $\\neq$ j, i $\\land$ j $\\in$ row, column or box, i $\\land$ j $\\notin$ known,\}\\\\ 0,\\text\{else,\}\\end\{cases\}\(10\)
with the following two conditions:

- •wu<0w\_\{u\}<0,
- •\|wu\|<1\|w\_\{u\}\|<1\.

While the weight matrix in Eq\.[10](https://arxiv.org/html/2607.15814#S3.E10)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](https://arxiv.org/html/2607.15814#S3.E8)to the following form:

ϕi˙=−funknown​\(ϕi\)​1N​∑jWi​j​sin⁡\(\(K2​\(ϕi−ϕj\)\)\)−funknown​\(ϕi\)​KS​sin⁡\(\(K​ϕi\)\),\\begin\{split\}\\dot\{\\phi\_\{i\}\}=&\-f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\\left\(\\frac\{K\}\{2\}\(\\phi\_\{i\}\-\\phi\_\{j\}\)\\right\)\}\\\\ &\-f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)K\_\{S\}\\sin\{\(K\\phi\_\{i\}\)\},\\end\{split\}\(11\)
with

funknown​\(ϕi\)=\{1,ϕi∈unknown,0,else\.f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)=\\begin\{cases\}1,\\text\{$\\phi\_\{i\}\\in$ unknown,\}\\\\ 0,\\text\{else\.\}\\end\{cases\}\(12\)

#### III\-A3Ensuring different phases

While the model presented in Eq\.[11](https://arxiv.org/html/2607.15814#S3.E11)accordingly describes the ONN\-ODE for Max\-K\-Cut, it has a serious shortcoming in regards to Sudoku\. Figure[2](https://arxiv.org/html/2607.15814#S3.F2)showcases an instance of a simple 4×\\times4 Sudoku puzzle with a violation of the rules\.

121123223Figure 2:Example of a 4×\\times4 Sudoku grid with a violation of the rules, where the thick numbers indicate the known values\.Clearly, the digit11in the upper left hand side corner violates two rules simultaneously, as there is already a known digit11in the same column as well as box\. However, despite breaking the rules, the ODE presented in Eq\.[11](https://arxiv.org/html/2607.15814#S3.E11)would produce a vanishing gradient, for the digits11and11are indeed2​π​nK\\frac\{2\\pi n\}\{K\}apart \(n=0n=0is also permitted\) and both values are one of the valid roots\. In summary, despite the structure, the ODE presented in Eq\.[11](https://arxiv.org/html/2607.15814#S3.E11)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 Max\-K\-Cut reads as follows:

ϕi˙=−funknown​\(ϕi\)​1N​∑jWi​j​sin⁡\(\(K2​\(ϕi−ϕj\)\)\)−funknown​\(ϕi\)​KS​sin⁡\(\(K​ϕi\)\)\+funknown​\(ϕi\)​KG​∑jWi​j′​g​\(ϕi,ϕj\),\\begin\{split\}\\dot\{\\phi\_\{i\}\}=&\-f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)\\frac\{1\}\{N\}\\sum\_\{j\}W\_\{ij\}\\sin\{\\left\(\\frac\{K\}\{2\}\(\\phi\_\{i\}\-\\phi\_\{j\}\)\\right\)\}\\\\ &\-f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)K\_\{S\}\\sin\{\(K\\phi\_\{i\}\)\}\\\\ &\+f\_\{\\text\{unknown\}\}\(\\phi\_\{i\}\)K\_\{G\}\\sum\_\{j\}W\_\{ij\}^\{\{\}^\{\\prime\}\}g\(\\phi\_\{i\},\\phi\_\{j\}\),\\end\{split\}\(13\)
withKGK\_\{G\}being a scaling factor,Wi​j′W\_\{ij\}^\{\{\}^\{\\prime\}\}a weight matrix that does not differentiate between knowns and unknowns, i\.e\.,

Wi​j′=\{−1,i≠j, i∧j∈row, column or box,0,else,W\_\{ij\}^\{\{\}^\{\\prime\}\}=\\begin\{cases\}\-1,\\text\{i $\\neq$ j, i $\\land$ j $\\in$ row, column or box,\}\\\\ 0,\\text\{else,\}\\end\{cases\}\(14\)
andg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)producing a non\-vanishing value if and only ifϕi\\phi\_\{i\}andϕj\\phi\_\{j\}are thennmultiples of2​π2\\pifar apart from each \(i\.e\. they output the same digit\) or in mathematical terms:

g​\(ϕi,ϕj\)=\{1,\|ϕi−ϕj\|≈2​π​n,0,else\.g\(\\phi\_\{i\},\\phi\_\{j\}\)=\\begin\{cases\}1,\|\\phi\_\{i\}\-\\phi\_\{j\}\|\\approx 2\\pi n,\\\\ 0,\\text\{else\.\}\\end\{cases\}\(15\)
We emulateg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)with the following function:

g​\(ϕi,ϕj\)=exp⁡\(\(−\(1−cos⁡\(\(ϕi−ϕj\)\)\)σ\)\),g\(\\phi\_\{i\},\\phi\_\{j\}\)=\\exp\{\\left\(\-\\frac\{\(1\-\\cos\{\(\\phi\_\{i\}\-\\phi\_\{j\}\)\}\)\}\{\\sigma\}\\right\)\},\(16\)
whereσ\\sigmatunes the width of the peaks around thenn\-th multiples of2​π2\\pi\.

### III\-BBenchmark set\-up

The model presented in Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13)has four tuneable parameters, namely:

- •the weight between two unknownswuw\_\{u\},
- •the scaling factor of theKK\-th harmonicsKSK\_\{S\},
- •the scaling factor of theg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)functionKGK\_\{G\},
- •the width of theg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)function’s peaksσ\\sigma\.

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 Py\-Sudoku library\[[36](https://arxiv.org/html/2607.15814#bib.bib9)\]\. In particular, we benchmark our proposed ONN\-Suduku solver \(see Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13)\) on Sudoku puzzles of sizes4×44\\times 4and9×99\\times 9for 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\]\[0\.125,0\.25,0\.375,0\.5,0\.625,0\.75,0\.875\]against an HNN Sudoku\-solver\[[27](https://arxiv.org/html/2607.15814#bib.bib10)\]and an already established ONN Sudoku\-solver\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]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 pre\-established roots\. All benchmarks are performed in Python\.

![Refer to caption](https://arxiv.org/html/2607.15814v1/x2.png)\(a\)
![Refer to caption](https://arxiv.org/html/2607.15814v1/x3.png)\(b\)

Figure 3:Accuracy as a function of unknown ratio of the HNN\-Sudoku solver\[[27](https://arxiv.org/html/2607.15814#bib.bib10)\], the conventional ONN\-Sudoku solver\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]and our proposed ONN\-Sudoku solver for \(a\)4×44\\times 4and \(b\)9×99\\times 9Sudoku 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\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]\.

## IVResults

First, we start with the accuracy benchmarks for4×44\\times 4and9×99\\times 9Sudoku puzzles\. In particular, we report on the best accuracies achieved with our ONN\- and the HNN\-Sudoku solver\[[27](https://arxiv.org/html/2607.15814#bib.bib10)\]for seven different unknown ratios\[0\.125,0\.25,0\.375,0\.5,0\.625,0\.75,0\.875\]\[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\]\[0\.1,0\.2,0\.3,0\.4,0\.5\]from the authors in\[[12](https://arxiv.org/html/2607.15814#bib.bib11)\]\. Figures[3\(a\)](https://arxiv.org/html/2607.15814#S3.F3.sf1)and[3\(b\)](https://arxiv.org/html/2607.15814#S3.F3.sf2)summarize all these accuracies as a function of the unknown ratio for sizes4×44\\times 4and9×99\\times 9, respectively\.

Let us start by analyzing the results for the size4×44\\times 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 than60%60\\%\. The accuracy drastically improves with the conventional ONN\-Sudoku solver achieving almost80%80\\%for an unknown ratio of50%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 of90%90\\%\. We speculate that this vastly superior performance can be traced back to the additional term presented in Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13)as it ensures a steady state if and only if all cells are correctly filled out\.

Let us now move on to the benchmarks on9×99\\times 9Sudoku 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 than40%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 of50%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 to25%25\\%and an accuracy of more than80%80\\%is obtained for an unknown ratio of37\.5%37\.5\\%, the precision starts to decrease, achieving approximately50%50\\%for an unknown ratio of50%50\\%and reaching the limits of the model at unknown ratios greater than70%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](https://arxiv.org/html/2607.15814#S3.E13)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 for4×44\\times 4and9×99\\times 9Sudoku puzzles\. In particular, we report on the order parameters associated with the best accuracies reached by our solver for seven different unknown ratios\[0\.125,0\.25,0\.375,0\.5,0\.625,0\.7,0\.875\]\[0\.125,0\.25,0\.375,0\.5,0\.625,0\.7,0\.875\]on 1000 Sudoku puzzles\. Figures[4\(a\)](https://arxiv.org/html/2607.15814#S4.F4.sf1)and[4\(b\)](https://arxiv.org/html/2607.15814#S4.F4.sf2)summarize all these order parameters as histograms for sizes4×44\\times 4and9×99\\times 9, respectively\.

![Refer to caption](https://arxiv.org/html/2607.15814v1/x4.png)\(a\)
![Refer to caption](https://arxiv.org/html/2607.15814v1/x5.png)\(b\)

Figure 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\]\[0\.125,0\.25,0\.375,0\.5,0\.625,0\.75,0\.875\]as well as \(a\)4×44\\times 4and \(b\)9×99\\times 9Sudoku puzzles, respectively\.Clearly, for4×44\\times 4Sudokus, 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, for9×99\\times 9Sudokus, 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 of0\.1250\.125, almost all order parameters are clustered around11, for an unknown ratio of0\.8750\.875, the vast majority of order parameters is below0\.50\.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 theg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)function in Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13), 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\.

## VDiscussion

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 for4×44\\times 4Sudokus and significantly improved accuracies for9×99\\times 9Sudokus\. However, it needs to be acknowledged that the model we considered \(see Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13)\) 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 for9×99\\times 9and unknown ratios larger than50%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 theg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)function \(see Eq\.[13](https://arxiv.org/html/2607.15814#S3.E13)\) is not scaled properly and/or the ONN did not have enough oscillatory cycles to settle\. Additionally, it may be worth looking into alternativeg​\(ϕi,ϕj\)g\(\\phi\_\{i\},\\phi\_\{j\}\)function formulations than the one presented in Eq\.[16](https://arxiv.org/html/2607.15814#S3.E16)\.

While our model deploys simple, sinusoidal Kuramoto\-ONNs\[[20](https://arxiv.org/html/2607.15814#bib.bib6),[43](https://arxiv.org/html/2607.15814#bib.bib5)\], there has recently been a push in the community to harness the dynamics of highly nonlinear oscillators for computing\[[33](https://arxiv.org/html/2607.15814#bib.bib12),[35](https://arxiv.org/html/2607.15814#bib.bib13),[25](https://arxiv.org/html/2607.15814#bib.bib14)\]\. 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\[[41](https://arxiv.org/html/2607.15814#bib.bib15),[8](https://arxiv.org/html/2607.15814#bib.bib16),[28](https://arxiv.org/html/2607.15814#bib.bib17),[16](https://arxiv.org/html/2607.15814#bib.bib18),[19](https://arxiv.org/html/2607.15814#bib.bib19),[29](https://arxiv.org/html/2607.15814#bib.bib20)\]\.

## VIConclusions

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 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 on4×44\\times 4and9×99\\times 9puzzles 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 on4×44\\times 4as well as9×99\\times 9puzzles, yielding almost flawless accuracies for4×44\\times 4Sudokus, while hitting a limit on9×99\\times 9Sudokus 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\.

## References

- \[1\]M\. Abernot, T\. Gil, M\. Jiménez, J\. Núñez, M\. J\. Avellido, B\. Linares\-Barranco, T\. Gonos, T\. Hardelin, and A\. Todri\-Sanial\(2021\)Digital implementation of oscillatory neural network for image recognition applications\.Frontiers in Neuroscience15,pp\. 713054\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[2\]M\. Abernot and A\. Todri\-Sanial\(2023\)Training energy\-based single\-layer hopfield and oscillatory networks with unsupervised and supervised algorithms for image classification\.Neural Computing and Applications35\(25\),pp\. 18505–18518\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[3\]P\. Bangert\(2012\)Optimization for industrial problems\.Springer Science & Business Media\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[4\]M\. K\. Bashar, Z\. Li, V\. Narayanan, and N\. Shukla\(2024\)An fpga\-based max\-k\-cut accelerator exploiting oscillator synchronization model\.IEEE\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[5\]Y\. Cheng, M\. Khairul Bashar, N\. Shukla, and Z\. Lin\(2024\)A control theoretic analysis of oscillator ising machines\.Chaos: An Interdisciplinary Journal of Nonlinear Science34\(7\)\.Cited by:[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.11)\.
- \[6\]D\. Efnusheva, A\. Cholakoska, and A\. Tentov\(2017\)A survey of different approaches for overcoming the processor\-memory bottleneck\.International Journal of Computer Science and Information Technology9\(2\),pp\. 151–163\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[7\]M\. Fathi, M\. Khakifirooz, and P\. M\. Pardalos\(2019\)Optimization in large scale problems: industry 4\.0 and society 5\.0 applications\.Vol\.152,Springer\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[8\]R\. FitzHugh\(1961\)Impulses and physiological states in theoretical models of nerve membrane\.Biophysical journal1\(6\),pp\. 445–466\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[9\]Y\. E\. Gonul, C\. E\. Kayan, I\. Mustafazade, N\. Kandasamy, and B\. Taskin\(2025\)GPU\-accelerated simulated oscillator ising/potts machine solving combinatorial optimization problems\.External Links:2505\.22631,[Link](https://arxiv.org/abs/2505.22631)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[10\]A\. Gower\(2025\)Learning at the speed of physics: equilibrium propagation on oscillator ising machines\.External Links:2510\.12934,[Link](https://arxiv.org/abs/2510.12934)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[11\]B\. F\. Haverkort and A\. Todri\-Sanial\(2026\)Overcoming quadratic hardware scaling for a fully connected digital oscillatory neural network\.Frontiers in NeuroscienceVolume 19 \- 2025\.External Links:[Link](https://www.frontiersin.org/journals/neuroscience/articles/10.3389/fnins.2025.1658490),[Document](https://dx.doi.org/10.3389/fnins.2025.1658490),ISSN 1662\-453XCited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[12\]B\. F\. Haverkort, F\. Sbravati, S\. Porfir, and A\. Todri\-Sanial\(2025\)Solving sudoku using oscillatory neural networks\.Neuromorphic Computing and Engineering\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§I](https://arxiv.org/html/2607.15814#S1.p4.7),[§I](https://arxiv.org/html/2607.15814#S1.p5.2),[Figure 3](https://arxiv.org/html/2607.15814#S3.F3),[Figure 3](https://arxiv.org/html/2607.15814#S3.F3.4.2),[§III\-B](https://arxiv.org/html/2607.15814#S3.SS2.p3.3),[§IV](https://arxiv.org/html/2607.15814#S4.p1.6)\.
- \[13\]J\. Hopfield\(1982\-05\)Neural networks and physical systems with emergent collective computational abilities\.Proceedings of the National Academy of Sciences of the United States of America79,pp\. 2554–8\.External Links:[Document](https://dx.doi.org/10.1073/pnas.79.8.2554)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[14\]F\. C\. Hoppensteadt and E\. M\. Izhikevich\(2000\)Pattern recognition via synchronization in phase\-locked loop neural networks\.IEEE Transactions on Neural Networks11\(3\),pp\. 734–738\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.3)\.
- \[15\]E\. Ising\(1924\)Beitrag zur theorie des ferro\-und paramagnetismus\.Ph\.D\. Thesis,Grefe & Tiedemann Hamburg, Germany\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3)\.
- \[16\]E\. M\. Izhikevich\(2003\)Simple model of spiking neurons\.IEEE Transactions on neural networks14\(6\),pp\. 1569–1572\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3),[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[17\]T\. Kanao and H\. Goto\(2022\-06\)Simulated bifurcation assisted by thermal fluctuation\.Communications Physics5\(1\),pp\. 153\(en\)\.External Links:ISSN 2399\-3650,[Link](https://www.nature.com/articles/s42005-022-00929-9),[Document](https://dx.doi.org/10.1038/s42005-022-00929-9)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[18\]R\. M\. Karp\(2009\)Reducibility among combinatorial problems\.In50 Years of Integer Programming 1958\-2008: from the Early Years to the State\-of\-the\-Art,pp\. 219–241\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1),[§III\-A1](https://arxiv.org/html/2607.15814#S3.SS1.SSS1.p1.1)\.
- \[19\]I\. Kovacic and M\. J\. Brennan\(2011\)The duffing equation: nonlinear oscillators and their behaviour\.John Wiley & Sons\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[20\]Y\. Kuramoto\(1984\)Chemical turbulence\.Springer\.Cited by:[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.11),[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[21\]A\. Lucas\(2014\)Ising formulations of many np problems\.Frontiers in physics2,pp\. 74887\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3),[§III\-A1](https://arxiv.org/html/2607.15814#S3.SS1.SSS1.p1.1)\.
- \[22\]O\. Maher, F\. Horst, F\. Sabo, D\. Fancheng, V\. Bragaglia, S\. Karg, M\. Sousa, A\. Todri\-Sanial, and B\. J\. Offrein\(2024\)Implementation of fitzhugh\-nagumo neurons using nanoscale vo 2 devices\.IEEE\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[23\]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\(2024\-04\)A CMOS\-compatible oscillation\-based VO2 Ising machine solver\.Nature Communications15\(1\),pp\. 3334\(en\)\.External Links:ISSN 2041\-1723,[Link](https://www.nature.com/articles/s41467-024-47642-5),[Document](https://dx.doi.org/10.1038/s41467-024-47642-5)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[24\]A\. Mallick, M\. K\. Bashar, Z\. Lin, and N\. Shukla\(2022\)Computational models based on synchronized oscillators for solving combinatorial optimization problems\.Physical Review Applied17\(6\),pp\. 064064\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§III\-A1](https://arxiv.org/html/2607.15814#S3.SS1.SSS1.p2.5),[§III\-A1](https://arxiv.org/html/2607.15814#S3.SS1.SSS1.p4.4)\.
- \[25\]O\. V\. Maslennikov, D\. S\. Shchapin, and V\. I\. Nekorkin\(2025\)Binary classification via spatiotemporal dynamics in reservoir computing rings of fitzhugh–nagumo neurons\.The European Physical Journal Special Topics234\(15\),pp\. 4033–4041\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[26\]T\. Miyato, S\. Löwe, A\. Geiger, and M\. Welling\(2025\)Artificial kuramoto oscillatory neurons\.External Links:2410\.13821,[Link](https://arxiv.org/abs/2410.13821)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[27\]V\. Mladenov, P\. Karampelas, C\. Pavlatos, and E\. Zirintsis\(2011\)Solving sudoku puzzles by using hopfield neural networks\.Proc\. of ICACM11,pp\. 174–179\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p4.7),[§I](https://arxiv.org/html/2607.15814#S1.p5.2),[Figure 3](https://arxiv.org/html/2607.15814#S3.F3),[Figure 3](https://arxiv.org/html/2607.15814#S3.F3.4.2),[§III\-B](https://arxiv.org/html/2607.15814#S3.SS2.p3.3),[§IV](https://arxiv.org/html/2607.15814#S4.p1.6)\.
- \[28\]J\. Nagumo, S\. Arimoto, and S\. Yoshizawa\(1962\)An active pulse transmission line simulating nerve axon\.Proceedings of the IRE50\(10\),pp\. 2061–2070\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[29\]E\. Panteley, A\. Loria, and A\. El Ati\(2015\)On the stability and robustness of stuart\-landau oscillators\.IFAC\-PapersOnLine48\(11\),pp\. 645–650\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[30\]H\. Ramsaueret al\.\(2021\)Hopfield networks is all you need\.External Links:2008\.02217Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[31\]N\. R\. Rohan, C\. Vigneswaran, S\. Ghosh, K\. Rajendran, A\. Gaurav, and V\. S\. Chakravarthy\(2025\)Deep oscillatory neural network\.Scientific Reports15\(1\),pp\. 40968\.External Links:[Document](https://dx.doi.org/10.1038/s41598-025-24837-4),[Link](https://doi.org/10.1038/s41598-025-24837-4)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[32\]F\. Rosenblatt\(1958\)The perceptron: a probabilistic model for information storage and organization in the brain\.\.Psychological review65\(6\),pp\. 386\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3)\.
- \[33\]F\. Sabo, F\. Du, N\. Dinç, and A\. Todri\-Sanial\(2026\)Harnessing dynamics of van der pol oscillators for phase computing with oscillatory neural networks\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[34\]F\. Sabo and A\. Todri\-Sanial\(2024\)ClassONN: classification with oscillatory neural networks using the kuramoto model\.IEEE\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[35\]M\. R\. E\. U\. Shougat, X\. Li, T\. Mollik, and E\. Perkins\(2021\)An information theoretic study of a duffing oscillator array reservoir computer\.Journal of Computational and Nonlinear Dynamics16\(8\),pp\. 081004\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[36\]J\. Sieu\(2024\)Py\-sudoku\.GitHub\.Note:\\urlhttps://github\.com/jeffsieu/py\-sudokuCited by:[§III\-B](https://arxiv.org/html/2607.15814#S3.SS2.p3.3)\.
- \[37\]L\. Sun, M\. X\. Burns, and M\. C\. Huang\(2025\-10\)General oscillator\-based Ising\-machine models with phase\-amplitude dynamics and polynomial interactions\.Physical Review Applied24\(4\),pp\. 044035\(en\)\.External Links:ISSN 2331\-7019,[Link](https://link.aps.org/doi/10.1103/7s4v-4hs4),[Document](https://dx.doi.org/10.1103/7s4v-4hs4)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[38\]A\. Todri\-Sanial, C\. Delacour, M\. Abernot, and F\. Sabo\(2024\)Computing with oscillators from theoretical underpinnings to applications and demonstrators\.Npj unconventional computing1\(1\),pp\. 14\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1),[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3),[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.11)\.
- \[39\]R\. Toews\(2020\-07\-17\)Deep learning’s carbon emissions problem\.Forbes\.External Links:[Link](https://www.forbes.com/sites/robtoews/2020/06/17/deep-learnings-climate-change-problem/)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[40\]S\. K\. Vadlamani, T\. P\. Xiao, and E\. Yablonovitch\(2024\-04\)Combinatorial optimization using the Lagrange primal\-dual dynamics of parametric oscillator networks\.Physical Review Applied21\(4\),pp\. 044042\(en\)\.External Links:ISSN 2331\-7019,[Link](https://link.aps.org/doi/10.1103/PhysRevApplied.21.044042),[Document](https://dx.doi.org/10.1103/PhysRevApplied.21.044042)Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1)\.
- \[41\]B\. van der Pol\(1926\)LXXXVIII\. on “relaxation\-oscillations”\.The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science2\(11\),pp\. 978–992\.Cited by:[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[42\]J\. von Neumann\(1993\)First draft of a report on the edvac\.IEEE Annals of the History of Computing15\(4\),pp\. 27–75\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[43\]T\. Wang and J\. Roychowdhury\(2019\)OIM: oscillator\-based ising machines for solving combinatorial optimisation problems\.InUnconventional Computation and Natural Computation: 18th International Conference, UCNC 2019, Tokyo, Japan, June 3–7, 2019, Proceedings 18,pp\. 232–256\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3),[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.11),[§II\-B](https://arxiv.org/html/2607.15814#S2.SS2.p1.3),[§V](https://arxiv.org/html/2607.15814#S5.p3.1)\.
- \[44\]X\. Yang and S\. Koziel\(2011\)Computational optimization and applications in engineering and industry\.Vol\.359,Springer Science & Business Media\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[45\]B\. Yegnanarayana\(2009\)Artificial neural networks\.PHI Learning Pvt\. Ltd\.\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3)\.
- \[46\]A\. R\. Yıldız\(2009\)An effective hybrid immune\-hill climbing optimization approach for solving design and manufacturing optimization problems in industry\.Journal of Materials Processing Technology209\(6\),pp\. 2773–2780\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p1.1)\.
- \[47\]J\. Zou, Y\. Han, and S\. So\(2009\)Overview of artificial neural networks\.Artificial neural networks: methods and applications,pp\. 14–22\.Cited by:[§I](https://arxiv.org/html/2607.15814#S1.p2.1),[§II\-A](https://arxiv.org/html/2607.15814#S2.SS1.p1.3)\.