# Lagrange oscillatory neural networks for constraint satisfaction and optimization

Corentin Delacour, Bram Haverkort, Filip Sabo, Nadine Azemard, Aida Todri-Sanial

## ► To cite this version:

Corentin Delacour, Bram Haverkort, Filip Sabo, Nadine Azemard, Aida Todri-Sanial. Lagrange oscillatory neural networks for constraint satisfaction and optimization. Neuromorphic Computing and Engineering, 2025, 5 (4), pp.044004. 10.1088/2634-4386/ae0eab. lirmm-05347029

HAL Id: lirmm-05347029

https://hal-lirmm.ccsd.cnrs.fr/lirmm-05347029v1

Submitted on 4 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

PAPER

# Lagrange oscillatory neural networks for constraint satisfaction and optimization

Corentin Delacour<fnref n="1" />,*, Bram Haverkort<fnref n="2" />,†, Filip Sabo<fnref n="2" />,†, Nadine Azemard<fnref n="1" /> and Aida Todri-Sanial<fnref n="1" />,<fnref n="2" />,‡

E-mail: delacour@ucsb.edu

**Keywords:** constrained optimization, Lagrange multipliers, oscillatory Ising machine, physics-based computing, Max-3-SAT

OPEN ACCESS

RECEIVED

12 May 2025

REVISED

26 September 2025

ACCEPTED FOR PUBLICATION

2 October 2025

PUBLISHED

10 October 2025

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

Abstract

Physics-inspired computing paradigms are receiving renewed attention to enhance efficiency in compute-intensive tasks such as artificial intelligence and optimization. Similar to Hopfield neural networks, oscillatory neural networks (ONNs) minimize an Ising energy function that embeds the solutions of hard combinatorial optimization problems. Despite their success in solving unconstrained optimization problems, Ising machines still face challenges with *constrained* problems as they can become trapped in infeasible local minima. In this paper, we introduce a Lagrange ONN (LagONN) designed to escape infeasible states based on the theory of Lagrange multipliers. Unlike existing oscillatory Ising machines, LagONN employs additional Lagrange oscillators to guide the system towards feasible states in an augmented energy landscape, settling only when constraints are met. Taking the maximum satisfiability problem with three literals as a use case (Max-3-SAT), we harness LagONN's constraint satisfaction mechanism to find optimal solutions for random SATlib instances with up to 200 variables and 860 clauses, which provides a deterministic alternative to simulated annealing for coupled oscillators. We benchmark LagONN with SAT solvers and further discuss the potential of Lagrange oscillators to address other constraints, such as phase copying, which is useful in oscillatory Ising machines with limited connectivity.

1. Introduction

Physics-inspired approaches, such as Ising machines, are being actively explored as potential efficient
hardware solutions for data-centric applications in artificial intelligence and optimization [2]. Ising
machines are closely related to the seminal works from Hopfield [3] and Tank [4] in the 1980s, who
introduced analog neural networks capable of naturally solving NP-hard combinatorial optimization
problems, such as the traveling salesman problem (TSP). Given a list of cities and their coordinates, TSP
seeks a minimum-distance tour that visits every city exactly once and returns to the starting point. Hopfield
and Tank's original approach maps cities and their positions in the tour to analog neurons, and intercity
distances to synaptic connections, such that the problem instance is 'physically wired'. The neurons then
evolve in parallel to minimize an Ising-like energy that encodes the tour length.

Despite its elegance, the system state can become trapped in local energy minima corresponding to infeasible tours, e.g. when two cities are visited simultaneously [5]. Constraints can theoretically be enforced with sufficiently large penalty coefficients [6], but this often slows convergence, as the dominance of constraint terms in the energy function hinders the search for optimal solutions. While software implementations can enforce hard constraints by generating only feasible samples, hardware Ising machines typically implement *soft* constraints embedded in the energy function. Consequently, current approaches to solving constrained optimization problems with Ising machines rely on tuning the penalty strength to balance feasibility and solution accuracy [7-9].

<footnote n="1">¹ Microelectronics Department, LIRMM, University of Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="2">² Electrical Engineering Department, Eindhoven Technical University, Eindhoven, The Netherlands</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
<figure><img src="image_4.png" /></figure>

With this energy function, what strategies, other than using impractically large penalty parameters, can we leverage for constraint satisfaction? Due to their continuous energy landscape, analog relaxations may provide techniques beyond those found in binary systems. In particular, they fit well with the formalism from Lagrange for constrained optimization on continuous functions [10, 11], which is a common principle to various physics-based solvers [12]. By relaxing the Ising energy to a continuous energy landscape, it is possible to apply Lagrange's theory and prevent the system state from settling to infeasible solutions, and rather force it to escape infeasible energy minima using additional variables—Lagrange multipliers $\lambda$ [13]. Conceptually, the dynamics consist of a competition between a gradient descent along the original variable direction (to minimize a cost function) and an ascent along the $\lambda$ direction (to enforce the constraints). Another interesting property of continuous variables is the parallel system evolution, contrasting with the typical sequential nature of simulated annealing [14], based on Gibbs sampling and only allowing sequential moves between connected spins.

In this paper, we revisit the concept of Lagrange multipliers applied to coupled phase oscillators, which we refer to as Lagrange oscillatory neural network (LagONN), as illustrated in figure 1. This platform is motivated by recent advances in coupled oscillator systems [15], which have demonstrated promising results in solving hard combinatorial optimization problems in the analog domain, including the Maximum cut problem [16–20], satisfiability (SAT) [21], and the Maximum independent set problem [22]. Dense implementations of coupled oscillators using CMOS technology are already available [19, 23–25], while emerging platforms such as spintronic devices [26, 27] and transition- or bistable-based devices [28–31] promise further scaling.

The primary objective of this work is to demonstrate constraint satisfaction with coupled phase oscillators without relying on quadratic energy penalties [6], thereby avoiding the trade-off between feasibility and solution accuracy common in constrained optimization problems. As an initial demonstration, we select a pure satisfaction problem: the Max-3-SAT problem, whose hardware acceleration has been extensively explored in literature [32–37]. The objective is to find a Boolean assignment of variable $x$ that maximizes the number of TRUE clauses composed of three literals $C_m = l_1^m \lor l_2^m \lor l_3^m$ in the formula:

$$f_B = C_1 \land C_2 \land \dots \land C_{M-1} \land C_M. \quad (1)$$
With $I_j^m \in \{x_1, \dots, x_N, \bar{x}_1, \dots, \bar{x}_N\}$. Assessing the SAT of $f_B$ is NP-complete [1], thus any problem in NP can be reduced to 3-SAT in theory.

For Ising machines, the Max-3-SAT problem can be approached in at least two ways. The first is an unconstrained optimization problem where the goal is to find the ground state of a static Ising energy combining all the clause terms (figure 1(b)). The second, which is the focus of this work, is to maximize constraint satisfaction via Lagrange multipliers that dynamically enforce the clauses [13], as shown in figure 1(c) for LagONN.

As the Lagrange formalism provides a deterministic mechanism to escape *infeasible* states, this translates into a ground state search for pure satisfaction problems like Max-3-SAT for which other classical solvers can get stuck at local minima (figure 1(c)). Consequently, LagONN constitutes a deterministic alternative to stochastic algorithms like simulated annealing [14] for ground state search. This is an interesting feature for oscillatory-based Ising machines, as harnessing noise and annealing can be challenging in practice.

The paper is organized as follows. After discussing prior work on physics-inspired solvers, we introduce a Lagrange oscillatory neural network (LagONN) that maps a 3-SAT formula to interconnected oscillatory subcircuits. The circuit dynamics seek an optimal saddle point in the energy landscape that satisfies all constraints. Next, we benchmark LagONN against SASAT [38], a simulated annealing algorithm, and the SAT solvers WalkSAT [39] and GNSAT-N [35] on Max-3-SAT instances from the SATlib library up to 200 variables and 860 clauses [40]. Finally, we discuss the broader utility of Lagrange oscillators by exploring their potential application to other constrained optimization problems, including copying phases in hardware with limited connectivity.

## 2. Background

### 2.1. Physics-inspired solvers for combinatorial optimization

At the heart of many physics-inspired solvers for combinatorial optimization lies the idea of mapping the problem's cost function to an energy, often set as the Ising Hamiltonian and expressed as:

$$H = - \sum_{i<j} J_{ij} S_i S_j - \sum_i h_i S_i \qquad (2)$$

where $S_i = \pm 1$ are the spins, $J_{ij}$ the interaction coefficients between spins, and $h_i$ is the external field applied to spin $S_i$. Finding the ground states of $H$ is in general NP-hard [6]. Updating a spin value $S_i$ causes an energy change $\Delta H_i$ that can drive the search in the energy landscape. While a deterministic search based only on energy minimization will likely get stuck at local minima, most state-of-the-art approaches allow some energy increase with nonzero probability to escape local minima. The standard choice for the system probability distribution is the Boltzmann law $\pi = \exp(-H/T)/Z$ that assigns high probabilities to low-energy states, where $T$ is a temperature parameter and $Z$ is the partition function. Using Markov chain Monte Carlo algorithms such as Gibbs sampling, one can then sample from the Boltzmann distribution by updating each spin sequentially, based on its local input $I_i = \sum_j J_{ij} S_j + h_i$ [41].

Since the goal of combinatorial optimization is to find low-energy states, the simulated annealing algorithm [14] is a natural choice for hardware implementation on noisy devices such as memristor arrays [42] or coupled probabilistic bits (p-bits), which inherently generate probabilistic spin updates [43]. Implementing Gibbs sampling with synchronous systems typically imposes a sequential spin update, although this is not a limitation for asynchronous systems [44] or sparse hardware with high parallelism [45]. In this paper, the system state components evolve in parallel owing to dynamics determined by differential equations. Our model is deterministic, although our numerical integration scheme introduces some errors, as discussed in the appendix B.

At first glance, a deterministic search may seem unsuitable for exponentially large spaces as the system state follows a determined trajectory depending on initialization, which may take exponential time before reaching an optimal point. Yet this limitation also applies to simulated annealing, where finite cooling schedules can trap variables in suboptimal states [46]. Other adiabatic approaches propose to slowly anneal the coupling amplitudes $|J_{ij}|$ in real-time, effectively shaping the energy landscape [47, 48], or to induce bifurcation phenomena during the adiabatic evolution of non-linear Hamiltonian systems, as in the simulated bifurcation algorithm [49, 50]. In our approach, the weight amplitude $|J_{ij}|$ is rather fixed, but its phase $\theta_{ij}$ varies in real-time, which can also be seen as an adaptive energy landscape.

Finally, many deterministic solvers rely on tailored differential equations [49, 51–53], where attractors in the phase space correspond to solutions, or deliberately introduce chaos to mimic Boltzmann sampling [54, 55]. A key consideration for building a corresponding physical machine is whether the system has bounded or unbounded variables. In [51, 56], the authors propose an analog model that provides an exponential
speed-up, but at the cost of unbounded variables and thus potentially exponential power. In practice, variable growth is capped by the finite power supply [57], which can lead to exponential computation times on hard instances [52]. In this work, we focus on dynamics realizable with coupled phase oscillators, ensuring bounded circuit power since variables are phases and the synaptic amplitudes $|J_{ij}|$ remain fixed. Consequently, we do not anticipate exponential speed-up, as confirmed by our time-to-solution (TTS) study on Max-3-SAT.

## 2.2. Computing with analog ONNs

This work focuses on ONNs operating in the phase domain: assuming identical oscillator frequencies, information is encoded in the phases $\phi_i$ relative to a reference oscillator. Under this assumption, and with symmetric synaptic weights $J_{ij} = J_{ji}$, coupled oscillators minimize an Ising-like energy and are similar to analog Hopfield neural networks [58–60]. With sinusoidal interactions, the network energy takes the form of a two-dimensional XY Ising Hamiltonian [61]:

$$E = -\sum_{i<j} J_{ij} \cos(\phi_i - \phi_j) - \sum_i h_i \cos(\phi_i) \quad (3)$$

which is an analog relaxation of the Ising model (equation (2)), since $E=H$ when the phases are binary (multiple of $\pi$). The connection between coupled oscillators and the Ising model can be illustrated with two coupled oscillators in the absence of external fields. For a coupling $J < 0$, the energy $E = -J\cos(\phi_1 - \phi_2)$ is minimized when the phase difference is $\pi$ (corresponding to opposite Ising spins), whereas for $J > 0$, the minimum occurs when the phases align (corresponding to identical spins). ONNs typically reach these energy minima via gradient descent, governed by:

$$\dot{\phi} = -\nabla_{\phi} E. \quad (4)$$

While the continuous nature of phases in $E$ has been exploited in efficient Max-Cut solvers [62, 63], most studies of coupled oscillators focus on the binary Ising model, often introducing harmonic signals to binarize phases, as proposed in [17]. In this work, we do not employ an additional binarization mechanism. As we will show, the fixed points naturally tend towards binary phases as the number of constraints (3-SAT clauses) increases.

In [21], the authors extend the energy $E$ to hypergraphs, incorporating high-order interactions between variables, whereas the classical Ising model (equation (2)) includes only pairwise (quadratic) spin interactions. This extension allows certain NP-hard problems, such as SAT, to be mapped directly to hardware without requiring quadratization, which adds auxiliary oscillators [64]. In this work, we focus on this approach, specifically considering simultaneous interactions among four oscillators $k, l, m, n$ whose energy can be written as:

$$E_{klmn} = J_{kl} \cos(\phi_k - \phi_l + \phi_m - \phi_n) = \mathrm{Re}(J_{kl} \exp[i(\phi_m - \phi_n)] \exp[i(\phi_k - \phi_l)]). \quad (5)$$

Examining equation (5), the additional third- and fourth-order interactions ($m$ and $n$) can be interpreted as a complex synapse connecting $k$ and $l$ with amplitude $J_{kl}$ and synaptic phase $\theta_{kl} = \phi_m - \phi_n$, which conveys the high order information from $m$ and $n$. This phase difference can also be interpreted as a delay $\tau_{kl} = (\phi_m - \phi_n)/\omega_0$ where $\omega_0$ is the oscillating frequency. With fourth-order interactions, the three variables of a 3-SAT clause plus an additional Lagrange variable could interact simultaneously. We note that implementing such high-order interaction is nontrivial and leave it as future work, though some conceptual ideas are discussed in section 3.3.

## 2.3. Constrained optimization with Lagrange multipliers

This work leverages the concept of Lagrange multipliers that are at the heart of many approaches to constrained optimization [10, 11, 65, 66]. For a cost function $f(x)$ with constraints $g(x)$, these methods typically define a Lagrange function:

$$L(x, \lambda) = f(x) + \lambda^T g(x) \quad (6)$$

where the constraints are satisfied when $g(x) = 0$, and $\lambda$ is a new variable called the Lagrange multiplier. Similar to the penalty method [10], the Lagrange approach weights the constraints and adds them to the cost function. The minimum of $L$ with respect to $x$ then provides a lower bound for the optimal constrained solution $f(x^*)$, since $\min_x L \le L(x^*, \lambda) = f(x^*)$. The purpose of the Lagrange multiplier $\lambda$ is to close this gap by finding an optimal $\lambda^*$ such that $\min_x L(x, \lambda^*) = f(x^*)$, or at least to reduce the gap as much as possible,
yielding the tightest lower bound on the optimal value [65]. Because L is concave in λ [67], the optimal λ* can be obtained by maximizing $\min_x L$:

$$L^* = \max_{\lambda} \min_{x} L(x, \lambda). \qquad (7)$$

For continuous functions f and g, this optimization problem can be interpreted as a Lagrange neural network, as defined in [68], where neurons perform gradient descent along x to minimize L, while Lagrange neurons perform gradient ascent in the λ-direction to reduce the gap (equation (7)) and reach an optimal solution f(x*) that satisfies the constraints. The corresponding dynamics are:

$$\left\{ \begin{array}{l} \tau \dot{x} = -\nabla_x L \\ \tau_\lambda \dot{\lambda} = +\nabla_\lambda L = g(x) \end{array} \right. \qquad (8)$$

where τ and τλ control the speeds of minimization and maximization, respectively. Reaching the global optimum f(x*) requires an ideal minimizer, and in non-convex settings, gradient descent can be trapped in local minima, yielding suboptimal solutions. Nevertheless, it is important to note that the search for a constrained solution cannot stop until the constraints are satisfied, i.e. g(x) = 0. Hence, for satisfaction problems such as Max-3-SAT, which we explore in this paper, the search continues until all clauses are TRUE. Another geometric interpretation of the dynamics in equation (8) is that they seek a saddle point of L, a minimum in the x-direction and a maximum in the λ-direction. In general, the energy landscape is non-convex and L* is not necessarily a saddle point [67]. For LagONN, however, we show that for satisfiable instances, L* corresponds to a saddle point that satisfies all constraints and can be reached via the competitive dynamics defined in equation (8).

# 3. Results

## 3.1. Mapping 3-SAT to ONNs

In this work, we build on the approach of Nagamatu and Yanaru [13] to solve the Max-3-SAT problem (equation (1)) using a Lagrange neural network with phase-based oscillatory neurons. Before introducing Lagrange multipliers, we now describe how to map a 3-SAT clause to coupled oscillators. First, we express each clause Cᵢ as an Ising energy Hᵢ that equals 0 if and only if Cᵢ is satisfied. For 3-SAT, there are four possible clause types, which we map to Ising energies as follows:

$$\begin{align*} C_1 &= X \lor Y \lor Z \rightarrow H_1 = 1 + S_X S_Y + S_X S_Z + S_Y S_Z - (S_X + S_Y + S_Z) - S_X S_Y S_Z \\ C_2 &= \bar{X} \lor Y \lor Z \rightarrow H_2 = 1 - S_X S_Y - S_X S_Z + S_Y S_Z - (-S_X + S_Y + S_Z) + S_X S_Y S_Z \\ C_3 &= \bar{X} \lor \bar{Y} \lor Z \rightarrow H_3 = 1 + S_X S_Y - S_X S_Z - S_Y S_Z - (-S_X - S_Y + S_Z) - S_X S_Y S_Z \\ C_4 &= \bar{X} \lor \bar{Y} \lor \bar{Z} \rightarrow H_4 = 1 + S_X S_Y + S_X S_Z + S_Y S_Z + (S_X + S_Y + S_Z) + S_X S_Y S_Z \end{align*} \qquad (9)$$

where spins Sⱼ = ±1 are the binary Ising variables with Sⱼ = +1 corresponding to TRUE. Writing the truth table for each clause, we find that Cᵢ is TRUE when Hᵢ = 0 and Cᵢ is FALSE when Hᵢ = 8. The next step is to map these Ising energies onto phase-based ONNs. We propose relaxing the binary Ising Hamiltonians to complex variables defined as:

$$\begin{align*} H_1 \rightarrow Z_1 &= 1 + e^{i(\phi_X - \phi_Y)} + e^{i(\phi_X - \phi_Z)} + e^{i(\phi_Z - \phi_Y)} - (e^{i\phi_X} + e^{i\phi_Y} + e^{i\phi_Z}) - e^{i(\phi_X - \phi_Y + \phi_Z)} \\ H_2 \rightarrow Z_2 &= 1 - e^{i(\phi_X - \phi_Y)} - e^{i(\phi_X - \phi_Z)} + e^{i(\phi_Z - \phi_Y)} - (-e^{i\phi_X} + e^{i\phi_Y} + e^{i\phi_Z}) + e^{i(\phi_X - \phi_Y + \phi_Z)} \\ H_3 \rightarrow Z_3 &= 1 + e^{i(\phi_X - \phi_Y)} - e^{i(\phi_X - \phi_Z)} - e^{i(\phi_Z - \phi_Y)} - (-e^{i\phi_X} - e^{i\phi_Y} + e^{i\phi_Z}) - e^{i(\phi_X - \phi_Y + \phi_Z)} \\ H_4 \rightarrow Z_4 &= 1 + e^{i(\phi_X - \phi_Y)} + e^{i(\phi_X - \phi_Z)} + e^{i(\phi_Z - \phi_Y)} + (e^{i\phi_X} + e^{i\phi_Y} + e^{i\phi_Z}) + e^{i(\phi_X - \phi_Y + \phi_Z)}. \end{align*} \qquad (10)$$

We deliberately introduce phase differences, e.g. $\phi_X - \phi_Y$, to define a complex relaxation Zᵢ of the conventional coupled-oscillators energy function (equation (3)). By construction, the complex relaxation Zᵢ equals the Hamiltonian Hᵢ for binary phases $\phi_j = \pi(1 - Sⱼ)/2$. However, non-binary phases can exist such that Zᵢ = 0, making the assignment of phases to Ising spins ambiguous. This issue is not caused by the high-order interaction terms but is inherent to the two-dimensional nature of the XY Ising model. A similar situation arises when solving the Max-cut problem with coupled oscillator systems without forcing binarization [17, 20], raising the question of how to round phases to spins.

In this work, we do not face this issue, as for more than a few clauses, the phase fixed points tend to be binary due to an overdetermined system of 2M equations and N phase variables, where M > N and M is the number of 3-SAT clauses. Leveraging this property, we design a coupled-oscillator module that minimizes $|Zᵢ|$, so that when $Zᵢ = 0$, the phases $\phi_X$, $\phi_Y$, and $\phi_Z$ directly yield the Boolean values S_X, S_Y, and S_Z satisfying clause Cᵢ.
<figure><img src="image_5.png" /></figure>

## 3.2. Enforcing constraints with an oscillatory-based Lagrange multiplier

Constraining the oscillators to satisfy a clause $C_i$ with energy $Z_i$ can be expressed via the Lagrange function $L_i(\phi, \lambda)$ as:

$$L_i(\phi, \lambda) = \lambda_R \mathrm{Re}[Z_i] + \lambda_I \mathrm{Im}[Z_i] \qquad (11)$$

where $\lambda_R$ and $\lambda_I$ are the Lagrange multipliers for the constraints, satisfied when $\nabla_{\lambda} L_i = (\mathrm{Re}[Z_i], \mathrm{Im}[Z_i]) = 0 \Leftrightarrow Z_i = 0$. Note that there are only two constraints and no cost function $f$, in contrast to the general Lagrange multiplier formulation (equation (6)). This defines a pure constraint satisfaction problem with constraint function $g = (\mathrm{Re}[Z_i], \mathrm{Im}[Z_i]) = 0$.

There are two possible interpretations for the Lagrange multipliers, which give rise to two distinct types of LagONNs:

1. $\lambda_R$ and $\lambda_I$ are *synaptic elements*. This raises the question of how to implement synapses $\lambda_R$ and $\lambda_I$ that evolve in real-time while having a limited range in a physical system.

2. $\lambda_R$ and $\lambda_I$ are *oscillatory variables*, i.e. phase oscillators. In this interpretation, only oscillating neurons and synapses with fixed amplitude are involved, potentially simplifying hardware implementation.

Focusing on the second interpretation, a specific choice for $\lambda$ simplifies the Lagrange function $L_i$ (equation (11)). We consider a Lagrange oscillator with the same frequency as the other neurons, phase $\phi_\lambda$, and unit amplitude (figure 2(b)). In the 2D plane, its corresponding unit vector $\vec{u}_\lambda$ has coordinates $(\cos \phi_\lambda, \sin \phi_\lambda)$, which we assign to multipliers $(\lambda_R, \lambda_I)$. Consequently, the Lagrange function reduces to the dot product between $\vec{u}_\lambda$ and the vector $\vec{Z}_i = (\mathrm{Re}[Z_i], \mathrm{Im}[Z_i])$ as:

$$L_i(\phi, \phi_\lambda) = \vec{u}_\lambda \cdot \vec{Z}_i = \cos \phi_\lambda \mathrm{Re}[Z_i] + \sin \phi_\lambda \mathrm{Im}[Z_i]. \qquad (12)$$
This yields a compact form for $L_i$, expressed here for the first type of clause (i = 1):

$$L_1(\phi, \phi_\lambda) = \cos \phi_\lambda \\ - \cos(\phi_X - \phi_\lambda) - \cos(\phi_Y - \phi_\lambda) - \cos(\phi_Z - \phi_\lambda) \\ + \cos(\phi_X - \phi_Y - \phi_\lambda) + \cos(\phi_X - \phi_Z - \phi_\lambda) + \cos(\phi_Z - \phi_Y - \phi_\lambda) \\ - \cos(\phi_X - \phi_Y + \phi_Z - \phi_\lambda). \qquad (13)$$

We identify the energy function as that of four sinusoidal coupled oscillators $\phi_X, \phi_Y, \phi_Z, \phi_\lambda$ with high-order interactions up to fourth order, as indicated by the last term and previously introduced in equation (5). Notably, the synaptic amplitudes appear as fixed binary weights $\pm 1$ for the cosine terms, ensuring that the LagONN variables remain bounded, since only phases vary and are bounded by definition.

Before deriving the corresponding circuit, we now analyze the properties of the energy landscape $L_i$. Introducing an oscillatory Lagrange multiplier generates a saddle point in the energy landscape where the constraint $Z_i = 0$ is satisfied, as formalized in the next theorem.

**Theorem 1.** Let $L_i(\phi, \phi_\lambda) = \vec{u}_\lambda \cdot \vec{Z}_i$ then:

1. $L_i$ has a at least one saddle point $L_i(\phi^*, \phi_\lambda^*) = 0$ such that $L_i(\phi^*, \phi_\lambda) \le L_i(\phi^*, \phi_\lambda^*) \le L_i(\phi, \phi_\lambda^*)$.

2. Such saddle point satisfies the constraint $Z_i(\phi^*) = 0$.

The theorem is proven in appendix D. Since LagONN's energy landscape $L_i(\phi, \phi_\lambda)$ has at least one saddle point satisfying the constraint, we set the phase dynamics for $\phi$ such that it minimizes $L_i(\phi, \phi_\lambda)$, while the Lagrange oscillator with phase $\phi_\lambda$ maximizes $L_i(\phi, \phi_\lambda)$. This aims to reach a saddle point where $Z_i(\phi^*) = 0$, as shown in figure 2(b). Accordingly, we propose the following phase dynamics for clause $C_i$:

$$\left\{ \begin{array}{l} \tau \dot{\phi}_j = -\nabla_{\phi_j} L_i(\phi, \phi_\lambda) = -\vec{u}_\lambda \cdot \partial \vec{Z}_i / \partial \phi_j \\ \\ \tau_\lambda \dot{\phi}_\lambda = +\nabla_{\phi_\lambda} L_i(\phi, \phi_\lambda) = \vec{u}_\lambda' \cdot \vec{Z}_i \end{array} \right. \qquad (14)$$

where $\tau$ and $\tau_\lambda$ are the time constants for the standard and Lagrange oscillators, and $\vec{u}_\lambda' = (-\sin \phi_\lambda, \cos \phi_\lambda)$. The time constants define the relative speed between gradient descent and ascent, which we take as $\tau = \tau_\lambda$ to achieve faster convergence in simulations (see appendix F). The dynamics can be interpreted as a competition between $\vec{Z}_i$ and $\vec{u}_\lambda'$, with the desired outcome being $Z_i = 0$ ($L_i$'s optimal saddle point) as a compromise between the two competing dynamics.

Figure 2 shows examples of ONN dynamics without and with a Lagrange oscillator. Without the Lagrange oscillator, the phases evolve to minimize Re[$\vec{Z}_i$] and settle to an undesired fixed point where $Z_i = -2$ (figure 2(a)). However, the desired state is $Z_i = 0$ (equation (10)), which would satisfy the clause $C_i$ for binary phases. Introducing a Lagrange oscillator as defined in equation (13), with competitive dynamics, constrains the phases to reach a saddle point satisfying $Z_i = 0$. Figure 2(b) shows an example of LagONN for the same initialization, highlighting the complex dynamics arising from this competition. Ultimately, the system settles to an optimal saddle point where $Z_i = 0$. Note that the Lagrange function (equation (12)) is not a Lyapunov function for the system, as it can increase over time (see appendix D).

### 3.3. Modular LagONN architecture for 3-SAT formula

We now propose a circuit implementation for a LagONN clause $C_i$ with energy $L_i$ and the competitive dynamics previously introduced (equation (14)). A network implementing a single clause as in equation (14) for $C_1 = X \lor Y \lor Z$ is shown in figure 3(a). The network includes a reference oscillator to measure phases and apply an external field to the Lagrange oscillator, so that in practice $\phi_j$ corresponds to $\phi_j - \phi_{REF}$. The main source of LagONN complexity is the synaptic array, which consists of delayed and weighted signals. In figure 3(a), a synapse $S_{ij}$ connecting oscillators i and j with weight $W_{ij}$ and phase $\theta_{ij}$ is represented as $S_{ij} = W_{ij}e^{i\theta_{ij}}$. This supposes having a mechanism to delay the synaptic input by $\theta_{ij}/\omega_0$ in real-time, which is the consequence of the fourth-order interaction terms in the LagONN function (equation (13)).

We believe this scheme is compatible with mixed-signal ONNs [20, 24], which could use digital synchronization circuits such as latches or counters to propagate the phase information from the third and fourth oscillators. Another approach would be to modulate the synaptic current amplitude between two oscillators in real-time with the phases of two others. Such circuitry would likely employ analog amplifiers, compatible with a variety of analog oscillators, including spintronic-based [26] or relaxation oscillators [20, 28, 69].
<figure><img src="image_6.png" /><figcaption>Figure 3. (a) LagONN network for solving clause C₁. Rectangles denote complex synapses in the form Jᵢⱼ exp (iθᵢⱼ) where θᵢⱼ is implemented with a delay θᵢⱼ/ω₀ in practice. The connections for θᵢⱼ are not shown. (b) Modular LagONN architecture to program any 3-SAT formula f<sub>B</sub> with N variables and M clauses. Each rectangle corresponds to a clause subcircuit and contains a Lagrange oscillator. The clause selection signal programs the clause subcircuit to one of the four possible clauses. N × 1 multiplexers route oscillator input and output signals to each clause according to the desired Boolean formula f<sub>B</sub>.</figcaption></figure>

We now introduce the circuit for a larger Boolean formula f<sub>B</sub> with N Boolean variables and M clauses
C<sub>m</sub> = l<sub>1</sub><sup>m</sup> ∨ l<sub>2</sub><sup>m</sup> ∨ l<sub>3</sub><sup>m</sup> defined as:

$$
f_B = C_1 \wedge C_2 \wedge \dots \wedge C_{M-1} \wedge C_M \qquad (15)
$$

with $l_i^m \in \{x_1, \dots, x_N, \bar{x}_1, \dots, \bar{x}_N\}$. The 3-SAT instance can be mapped to LagONN modules where a module *m* corresponds to a clause *C*<sub>*m*</sub>, and each literal *l*<sub>*i*</sub><sup>*m*</sup> corresponds to an input port of that module (figure 3(b)). Each module implements the network shown in figure 3(a) and is highlighted by the blue box. Implementing the AND operation ‘**∧**’ between two clauses is unnecessary, as each Lagrange oscillator evolves to satisfy its respective clause. Consequently, the entire network maximizes the number of TRUE clauses in *f*<sub>B</sub>.

LagONN modules are connected as follows. If two clauses *m* and *n* share a literal at positions *k* and *l*: either identical $l_k^m = l_l^n$ or negated $l_k^m = \overline{l_l^n}$ corresponding to variable $\phi_{x_j}, j \in \{1, \dots, N\}$, we connect input ports *k* and *l* so that the synapses from both clauses influence $\phi_{x_j}$. Repeating this procedure for every pair of clauses ensures that each variable $\phi_{x_j}$ is influenced by all corresponding clauses, yielding the total LagONN function:

$$
L_T(\phi_x, \phi_\lambda) = \sum_{m=1}^{M} u_\lambda^m \cdot \vec{Z}_m(\phi_x) \quad (16)
$$

where $\phi_x = (\phi_{x_1}, \phi_{x_2}, \dots, \phi_{x_{N-1}}, \phi_{x_N})^T$ is the vector of phases corresponding to the Boolean variables and
$\phi_\lambda = (\phi_{\lambda_1}, \phi_{\lambda_2}, \dots, \phi_{\lambda_{M-1}}, \phi_{\lambda_M})^T$ contains all the Lagrange variables.

The proposed modular architecture offers two main advantages. First, it avoids the need for large
synaptic arrays, as synapses are confined within each module. Second, it preserves high-order synaptic
interactions locally within the modules. In contrast, a programmable fully-connected design with N + M
oscillators would require a (N + M)² synaptic array to support any f<sub>B</sub> instance. Propagating high-order
interactions across such a large array is particularly challenging, as it is not straightforward with a standard
two-dimensional grid layout. However, the modular approach introduces a different type of
two-dimensional array: since each module input/output can be connected to N different input/output lines,
a programmable architecture would require 6M N × 1 multiplexers. Thus, compared to a fully connected
design, there is a trade-off between the overhead from multiplexers and the complexity of propagating
high-order interactions. The comparison between the modular and fully connected LagONN architectures
scaling is summarized in table 1.
<table><caption>**Table 1.** Comparison between two possible LagONN architectures.</caption><thead><tr><td>LagONN architecture</td><td>Fully-connected</td><td>Modular</td></tr></thead><tbody><tr><td>Oscillators</td><td>N + M</td><td>N + M</td></tr><tr><td>Size of synaptic array</td><td>(N + M)²</td><td>4²/module (clause)</td></tr><tr><td>Distance of high-order interaction</td><td>N + M (global)</td><td>4 (local)</td></tr><tr><td>N × 1 multiplexers</td><td>0</td><td>6M</td></tr></tbody></table>

**ONN phase dynamics**

<figure><img src="image_7.png" /></figure>

**LagONN phase dynamics**

<figure><img src="image_8.png" /></figure>

**Cost evolution**

<figure><img src="image_9.png" /></figure>

**LagONN energy for first 4 clauses**

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Phase dynamics for the satisfiable SATlib instance 'cnf-20-01' with N = 20 variables and M = 91 clauses. (a) Standard ONN dynamics quickly converge toward a sub-optimal solution where 89/91 of clauses are satisfied. (b) With the same initialization, the Lagrange version takes more time to reach a fixed point. Reading out the phases gives an optimal Boolean assignment where all clauses are satisfied. (c) Cost function comparison between the two approaches. LagONN finds an assignment of optimal phase around the same time the standard ONN settles (≈10 oscillation cycles). By measuring the cost in real time, we can stop the run when a target cost is reached without waiting for convergence. (d) Dynamics for four LagONN energy terms Z_m corresponding to the first four clauses. While the dynamics almost seem chaotic, they evolve to reach a target saddle point where all Z_m = 0. Ultimately, the dynamics converge toward a fixed point at t = 300 oscillation cycles where all Z_m = 0.</figcaption></figure>

### 3.4. LagONN competitive dynamics

We now express the phase dynamics of the entire LagONN circuit shown in figure 3(b). Assuming that the formula $f_B$ is satisfiable, we extend theorem 1 to the total Lagrange function $L_T$ with M clauses.

**Theorem 2.** Let $L_T(\phi, \phi_\lambda) = \sum_m^M u_\lambda^m \vec{Z}_m$ and the formula $f_B = C_1 \bigwedge C_2 \bigwedge \dots \bigwedge C_{M-1} \bigwedge C_M$ is satisfiable, then:

1. $L_T$ has at least one saddle point $L_T(\phi^*, \phi_\lambda^*) = 0$ such that $L_T(\phi^*, \phi_\lambda) \le L_T(\phi^*, \phi_\lambda^*) \le L_T(\phi, \phi_\lambda^*)$.

2. Such saddle point satisfies the constraints $Z_m(\phi^*) = 0$ for all clauses.

The proof is in the appendix D. As for a single clause, we set the dynamics of the entire LagONN system to reach a saddle point as:

$$ \begin{cases} \tau \dot{\phi}_x = -\nabla_{\phi_x} L_T = -\sum_m u_\lambda^m \partial \vec{Z}_m / \partial \phi_x \\ \tau_\lambda \dot{\phi}_\lambda = +\nabla_{\phi_\lambda} L_T = \sum_m u_\lambda^m \partial \vec{Z}_m. \end{cases} \qquad (17) $$

Figures 4(a) and (b) show examples of phase dynamics without and with Lagrange oscillators for N = 20 variables and M = 91 clauses (satisfiable instance cnf-20-01 from SATlib [40]). Without Lagrange oscillators, the ONN phases converge to non-binary values. To extract Boolean assignments, we round each phase to the
nearest multiple of $\pi$, which yields a sub-optimal solution with two unsatisfied clauses. With Lagrange oscillators, the dynamics become more complex and take longer to settle. However, phases settle to multiples of $\pi$, eliminating the need for rounding, and providing an optimal Boolean assignment where all 91 clauses are satisfied. Since each clause enforces two constraints, $Re[Z] = 0$, and $Im[Z] = 0$, the system is overdetermined with $2M$ equations and only $N$ unknowns ($M > N$). We hypothesize that this high level of frustration restricts the saddle points to binary phases only.

Figure 4(d) shows the dynamics of the first four energy terms $Z_m$, which are simultaneously attracted to the origin $Z_m = 0$ corresponding to a target saddle point. An optimal assignment is obtained when all $Z_m$ trajectories converge to zero. Figure 4(c) compares the cost function evolution for the two cases. Interestingly, LagONN dynamics reduce the cost as quickly as the ONN, reaching about two unsatisfied clauses after a single oscillation cycle. However, unlike the ONN that gets trapped, LagONN continues to evolve and eventually finds an optimal Boolean assignment. Comparing the phase dynamics from figure 4(b) with the cost evolution in figure 4(c) reveals that some phases keep evolving with little effect on the cost. For this reason, in all the simulations reported in this paper, we monitor the cost in real-time and stop the simulation once the cost reaches a satisfactory value (set to 0 for satisfiable instances), as described in the appendix A. In other words, we do not require full convergence to the saddle point shown in figure 4(b), since its stability under the dynamics of equation (17) is not guaranteed (see appendix E).

### 3.5. Comparison with simulated annealing and TTS

We compare LagONN's search with a simulated annealing algorithm for SAT [38] detailed in the appendix C. Each SAT variable is sequentially flipped and the corresponding cost change $\delta$ is calculated. A flip is accepted with probability $1/(1 + \exp(\delta/T))$ where $T$ is a temperature parameter. The temperature decreases exponentially from $T_{max} = 1$ to $T_{min} = 0.01$, values chosen empirically by inspecting cost trajectories across different instance sizes. Following the exponential schedule proposed in [38], we obtain satisfactory results, though we do not claim optimality.

An example of simulated annealing run is shown in figure 5(a) (left) for 100 variables and 430 clauses, with the x-axis indicating the number of steps (tentative flips). The 100 variables are arbitrarily initialized to '1', and after the first 100 steps (one sweep), each variable has been considered once and potentially flipped according to the sigmoid probability at $T = T_{max}$. This initial sweep induces a sharp cost drop of about 30 satisfied clauses. As the temperature decreases, the algorithm enters a regime with fewer flips and the cost gradually declines. The run terminates once an optimal variable assignment is found.

A LagONN simulation is shown in figure 5(a) (right) for the same variable initialization and a random Lagrange oscillator initialization. The cost trajectory differs qualitatively: unlike simulated annealing, LagONN has no annealing schedule and continues making uphill moves even at low cost, which simulated annealing would rarely accept at low temperature. Crucially, LagONN dynamics do not halt until all constraints are satisfied, whereas simulated annealing can freeze into a suboptimal state within finite annealing time. However, this does not imply greater efficiency, as shown in the next TTS comparison.

We now compare LagONN and simulated annealing on hard random instances from SATlib [40] with increasing sizes $(N, M) \in \{(20, 91), (50, 218), (75, 325), (100, 430), (125, 538), (150, 645), (200, 860)\}$. These instances are part of a well-known benchmark, lying near the computational phase transition where the clause-to-variable ratio is $M/N \approx 4.3$ [70], making them especially challenging. For each $(N, M)$, we evaluate the first 100 instances for simulated annealing and the first 30 for LagONN (fewer instances due to the longer runtime of LagONN simulations). Each instance is run with 100 random initializations, where the phases are drawn uniformly, and simulations are executed for a fixed runtime $t_{max}$. From these runs, we estimate the success probability $p_s$. As in other Ising-machine benchmarks, we report the TTS, defined as the expected time to reach an optimal solution with probability 0.99. For each instance, the TTS is given by:

$$ \text{TTS} = t_{\text{max}} \frac{\log(0.01)}{\log(1 - p_s)}. \qquad (18) $$

We apply the same TTS formula for simulated annealing, where $t_{max} = n_{max}$ corresponds to the number of algorithmic steps per run (maximum number of variable updates). When analyzing TTS scaling with system size, it is crucial to optimize the runtime for each size. Otherwise, misleading effects may appear such as artificially flat curves when $t_{max}$ is set too large, or TTS overestimation when $t_{max}$ is too small for larger instances [71]. To mitigate these issues, we iteratively increase $t_{max}$ and $n_{max}$ with the system size, up to $t_{max} \le 10^6$ oscillation cycles, and $n_{max} \le 2 \times 10^7$ steps. This ensures that the success probability remains within the range $0.1 < p_s < 0.9$ for all tested instances.

Figure 5(b) shows the TTS scaling with the number of variables at a fixed clauses-to-variables ratio $M/N \approx 4.3$, plotted on a semilog scale with blue boxes corresponding to 1st, 2nd, and 3rd quartiles. Qualitatively, simulated annealing exhibits more favorable scaling than LagONN on this benchmark, as its
<figure><img src="image_11.png" /><figcaption>**Figure 5.** (a) Cost evolution comparison for a SATlib instance with 100 variables and 430 clauses, where all variables are initialized to 1. The temperature in simulated annealing decreases exponentially with the number of steps as described in the appendix C. (b) Estimation of the time to reach an optimal solution (TTS) with 99% probability for SATlib instances with a clauses-to-variables ratio $M/N \approx 4.3$. Boxes show the 1st, 2nd, and 3rd quartiles computed for the first 30 satisfiable instances from SATlib for LagONN, and 100 instances for simulated annealing. The error bars show the min-max values, and red circles are the averages. For both methods, we fit the logarithm of data to estimate the mean TTS scaling as $\sim \exp(aN + b)$ with $a = 0.059$, $b = 3.58$ for LagONN and $a = 0.055$, $b = 8.16$ for simulated annealing (red dashed lines). Simulated annealing has an advantageous scaling for the median TTS fitted as $\sim \exp(c\sqrt{N} + d)$ with $c = 0.719$, $d = 5.49$, whereas we fit LagONN's median TTS as $\sim \exp(aN + b)$ with $a = 0.056$ and $b = 2.43$ (blue dashed curves). (c) For a fixed number of variables $N = 50$, we vary the number of clauses $M$ from 100 to 218 and compute TTS for 100 instances per point. The median TTS scales exponentially with the number of clauses for both methods. (d) Benchmark of median TTS with SAT solvers from Pedretti et al [35] based on stochastic local search for satisfiable SATlib instances.</figcaption></figure>

slope flattens at larger sizes. Specifically, we fit the median TTS as $\sim \exp(c\sqrt{N} + d)$ with $c = 0.719$, $d = 5.49$, whereas LagONN's median TTS follows a straight line in the semilog plot, fitted as $\sim \exp(aN + b)$ with $a = 0.056$ and $b = 2.43$ (blue dashed curves). For simulated annealing, we excluded a polynomial scaling for the median TTS, as the log-log plot did not yield a straight line.

At large sizes, LagONN's TTS has fewer outliers than simulated annealing, with a mean TTS scaling similarly to the median as $\sim \exp(aN + b)$ with $a = 0.059$, $b = 3.58$. In contrast, simulated annealing sometimes encounters very hard instances, leading to more than five decades of variations at $N = 200$, which pushes the mean scaling to $\sim \exp(aN + b)$ with $a = 0.055$, $b = 8.16$, comparable in slope to LagONN. Overall, while simulated annealing achieves a more favorable median scaling, LagONN remains competitive due to its lower prefactor $\exp(2.43) < \exp(5.49)$, which keeps it competitive up to $N \le 200$.

As shown in figure 5(c) for $N = 50$, both algorithms exhibit an exponential increase in TTS as the clauses-to-variables ratio approaches the computational phase transition at $M/N \approx 4.3$. Between $M/N = 2$ and $4.3$, the TTS increases by roughly two orders of magnitude for simulated annealing and by about one order of magnitude for LagONN, confirming that problem difficulty grows with the number of clauses. Fixing $M/N = 2$ for all system sizes halves LagONN's median scaling slope in figure 5(b), yielding a new slope of $a = 0.027$ and a median TTS of $4.4 \times 10^3$ oscillations for $N = 200$. We remain cautious about these scaling estimates, as they are based on a limited number of instances and sizes, constrained by the very long simulation times at $M/N \approx 4.3$ (see appendix B). Running the same TTS scaling experiment at $M/N \approx 4.3$ with a reduced integration time step (linearly annealed from 0.15 to 0.015) did not significantly improve LagONN's exponential scaling.

### 3.6. Benchmark with SAT solvers

We benchmark LagONN with SAT solvers based on stochastic local search (SLS), namely WalkSAT [39, 72] and a hardware-enhanced version of GSAT, GNSAT-N [35]. SLS algorithms flip variables iteratively, similar to simulated annealing, but differ in how they select variables. One of the simplest methods is GSAT [73], a greedy algorithm that flips a variable which maximizes the overall cost reduction or gain. Pedretti et al [35]
proposes GNSAT-N, a GSAT hardware acceleration and improvement by adding normal-distributed noise to the gain, thereby allowing random walks. WalkSAT has another strategy to select variables and focuses on unsatisfied clauses. It first selects an unsatisfied clause and computes the number of new unsatisfied clauses (called breaks) induced when flipping each variable of the clause. If for some variable of the clause, break = 0 (zero-damage flip [72]), it is flipped. Otherwise, the best variable (minimizing breaks) is flipped with some probability, or a random variable is flipped (random walk). It is then a compromise between greedy moves and a random walk, set by the probability parameter.

Figure 5(d) shows the median TTS (steps) for WalkSAT and GNSAT-N reported by Pedretti et al [35] on uniform satisfiable SATlib instances, alongside LagONN's TTS expressed in number of oscillations. WalkSAT and GNSAT-N have similar TTS and exhibit a clear advantage in scaling over simulated annealing and LagONN. We interpret the performance gap between SAT solvers and simulated annealing as follows. Although the probability of flipping a variable in simulated annealing depends on the cost reduction (following a sigmoidal function), its current implementation (SASAT [38]) evaluates each variable iteratively in a fixed order. This contrasts with GSAT or WalkSAT, which have more sophisticated mechanisms to select variables or clauses before making a flip, potentially making more efficient moves.

Although LagONN also has a greedy mechanism to follow steepest trajectories minimizing the cost (gradient descent), its dynamics are strongly perturbed by the competitive forces introduced by the Lagrange oscillators, as shown in figure 2. This interplay of simultaneous gradient descent and ascent appears less efficient than the highly tuned SLS searches designed for SAT, which is reflected in LagONN's less favorable TTS scaling. Nevertheless, its lower prefactor keeps LagONN competitive for $N \le 150$, and the results overall confirm its ability to enforce constraint satisfaction by escaping infeasible states in a deterministic manner. More generally, these findings suggest that similar competitive oscillator dynamics could be applied to constrained optimization problems formulated through a Lagrange function (equation (6)) as we discuss next with the example of phase copying. In this case, Lagrange oscillators are used to enforce the constraints during optimization rather than serving purely as a search mechanism, as in Max-3-SAT.

## 4. Summary and discussion

In this work, we demonstrated how introducing additional Lagrange oscillators into coupled oscillator systems enables constraint satisfaction, exemplified by the Max-3-SAT problem. These problems represent special cases where the Lagrange function reduces to the constraint function $g$, which counts unsatisfied clauses, without any additional cost term $f$. Because the Lagrange oscillators drive the system until all constraints are satisfied, LagONN continues its search until it finds an optimal solution, a distinct advantage over simulated annealing, which can become trapped in local minima. Another benefit is that LagONN is noiseless and requires no annealing schedule, a practical advantage given that tuning noise in physical oscillator implementations can be challenging. On the other hand, LagONN's search for optimal Boolean assignments is less effective than simulated annealing and state-of-the-art SAT solvers, as indicated by our TTS analysis. Furthermore, the stability of the saddle points reached by the dynamics in equation (17) is not guaranteed, although stabilization mechanisms are possible (see appendix E). A side-by-side comparison of the strengths and limitations of each method is provided in table 2.

Beyond pure constraint satisfaction problems, LagONN can also be applied to more general constrained problems involving both a cost $f$ and constraints $g$. For instance, figure 6 illustrates the problem of phase copying, which can arise in large-scale hardware implementations. When a problem is represented by a dense graph that cannot be directly implemented due to the quadratic number of edges, one can introduce copy nodes that copy each other to distribute the edges [74]. This scenario is typical for quantum annealers with limited connectivity [75]. Here, the constraint $g$ enforces equality between copies, while the cost function $f$ corresponds to the Ising energy of the original graph. To copy the spin values, one can introduce a ferromagnetic coupling $J_c$ between nodes, corresponding to the penalty term $\frac{1}{2}(S_1 - S_2)^2 = -J_c S_1 S_2$. In principle, this positive coupling should be large enough to satisfy the constraint. However, in practice, too strong couplings can rigidify the dynamics and hinder the search for ground states [7, 74].

Figure 6(a) shows an example of a coupled oscillator network with standard coupling $J = -1$ (blue edges) and copy coupling $J_c$ (dashed black edges). The cost function $f$ to minimize is the XY Ising energy $E$ for phase oscillators (equation (3)), determined by the symmetric couplings. The copy constraints are defined as $g = (\cos\phi_1 - \cos\phi_2, \cos\phi_3 - \cos\phi_4) = 0$. When the constraints are satisfied, the graph reduces to a 4-node fully-connected graph. Figure 6(b) shows an example of phase dynamics with a weak copy coupling $J_c = +0.5$ where the constraints are not fully enforced, e.g. $\cos\phi_1 \neq \cos\phi_2$.

Instead of increasing $J_c$ at the risk of rigidifying the system, one can introduce a Lagrange oscillator for each constraint to enforce phase equality, as illustrated in figure 6(c). For example, for the phases $\phi_1$ and $\phi_2$, the corresponding copy constraint can be written as a vector $Z_{1-2} = \exp(i\phi_1) - \exp(i\phi_2) = 0$. Analogous to
<table><caption>**Table 2.** Comparison between simulated annealing and Lagrange ONN for solving the Max-3-SAT constraint satisfaction problem, from both conceptual and practical (hardware) points of view.</caption><thead><tr><th></th><th>Simulated annealing</th><th>Lagrange ONN</th></tr></thead><tbody><tr><td>Strengths</td><td>Advantageous TTS scaling.<br>Convergence at cold temperature.</td><td>Noiseless approach. Cannot settle into<br>infeasible states.</td></tr><tr><td>Weaknesses</td><td>Requires careful noise tuning. Can be<br>stuck at infeasible states.</td><td>Seemingly worse TTS scaling. Stability<br>is not guaranteed.</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 6.** LagONN for phase copying on an initial graph with negative coupling $J = -1$ (blue edge). (a) In practice at a large scale, Ising machines often have limited hardware connectivity, and dense graphs are mapped to a sparser network by introducing copy nodes [74]. Ideally, the copy nodes should copy the original values, e.g. $\cos \phi_2 = \cos \phi_1$. This is generally achieved with a ferromagnetic coupling between the two nodes ($J_c = +0.5$ here). In theory, the positive coupling should be large enough to keep the two variables identical. However, a large coupling value introduces rigidity in the system, hindering the ground state search. (b) For weak coupling values such as in this example, the copied phases are not equal ($\cos \phi_1 \neq \cos \phi_2$). (c) Lagrange oscillators can enforce the copy constraints between two copied nodes, as shown in this example for $J_{\lambda\phi} = 1$.</figcaption></figure>

the Lagrange function defined for 3-SAT (equation (16)), the Lagrange term for this constraint is expressed as:

$$ \vec{u}_{\lambda_1} \cdot \vec{Z}_{1-2} = \cos(\phi_1 - \phi_{\lambda_1}) - \cos(\phi_2 - \phi_{\lambda_1}) \quad (19) $$

where $\phi_{\lambda_1}$ is the phase of the Lagrange oscillator connected to the first pair of oscillators. The total Lagrange function for this example is then:

$$ L(\phi, \phi_{\lambda}) = E(\phi) + J_{\lambda\phi} \left( \vec{u}_{\lambda_1} \cdot \vec{Z}_{1-2} + \vec{u}_{\lambda_2} \cdot \vec{Z}_{3-4} \right) \quad (20) $$

where $J_{\lambda\phi}$ sets the constraint strength or equivalently, the speed of Lagrange oscillators. Implementing competitive dynamics seeking a saddle point in the Lagrange landscape as in equation (17), $\phi_{\lambda}$ performs gradient ascent until the constraints are satisfied, while other phases $\phi$ follow gradient descent to minimize $L$. The opposite gradient signs induce asymmetric couplings between oscillators and Lagrange oscillators, illustrated by the arrows in figure 6(a) with weights set to $J_{\lambda\phi} = -J_{\phi\lambda} = 1$.

A key limitation of this approach is that the stability of fixed points is not guaranteed and requires further analysis. For example, in simulations, increasing the strength $J_{\lambda\phi}$ speeds up the Lagrange oscillators and introduces transient and decaying phase oscillations. In contrast, increasing the copying strength $J_c$ acts like additional damping and smooths the dynamics. Combining energy penalties and Lagrange multipliers could enhance existing oscillatory Ising machines, as proposed in previous work for knapsack problems [76], potentially yielding orders-of-magnitude speed-ups. A detailed exploration of these regimes is left for future work.

## 5. Conclusion

This article introduced LagONN, a LagONN that enforces constraint satisfaction as demonstrated for the Max-3-SAT problem. Unlike gradient-descent-based approaches such as oscillatory Ising machines, which can be trapped in infeasible local minima, LagONN employs additional Lagrange oscillators to ensure that 3-SAT clauses are satisfied. Conceptually, these Lagrange variables provide alternative pathways in the energy landscape to escape local minima and reach optimal states where phases correspond to optimal Boolean values. When benchmarked on SATlib instances up to 200 variables and 860 clauses against simulated annealing and SAT solvers, LagONN simulations exhibited less favorable TTS scaling with increasing
problem size, yet remained competitive for the tested problem sizes ($N \le 200$) due to a smaller prefactor in the scaling function. For the Max-3-SAT constraint satisfaction problem, LagONN offers a deterministic search method that eliminates the need for careful noise control required in simulated annealing, which can be difficult to implement in hardware. Furthermore, the example of phase copying illustrates how oscillatory-based Ising machines augmented with Lagrange oscillators can enforce constraint satisfaction in general optimization tasks, going beyond traditional penalty methods.

## Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/corentindelacour/Lagrange-oscillatory-neural-network.

## Acknowledgments

This work was supported by the European Union's Horizon 2020 research and innovation program, EU H2020 NEURONN (www.neuronn.eu) project under Grant 871501. B M, F S and A T S acknowledge support from the European Union's Horizon Europe research and innovation programme, PHASTRAC project under Grant Agreement No. 101092096 and Dutch Research Council's AiNed Fellowship research programme, AI-on-ONN project under Grant Agreement No. NGF.1607.22.016, as well as funding from the European Research Council ERC THERMODON project under Grant Agreement No. 101125031.

## Code availability

Matlab codes and data are available at the following GitHub repository: https://github.com/corentindelacour/Lagrange-oscillatory-neural-network

## Author contributions

Coretin Delacour ID 0000-0001-6386-9588

Conceptualization (lead), Formal analysis (lead), Software (lead), Writing – original draft (lead)

Bram Haverkort ID 0009-0009-4391-229X

Data curation (supporting), Investigation (supporting), Writing – original draft (supporting)

Filip Sabo ID 0009-0003-7885-9747

Data curation (supporting), Investigation (supporting), Writing – original draft (supporting)

Nadine Azemard

Validation (supporting), Writing – review & editing (supporting)

Aida Todri-Sanial ID 0000-0001-8573-2910

Conceptualization (supporting), Funding acquisition (lead), Project administration (lead), Resources (lead), Supervision (supporting), Writing – review & editing (supporting)

## Appendix A. LagONN cost function

Here, we describe how LagONN's cost is monitored during the simulation. From the SATlib.cnf files [40], we build a system of differential equations for each instance (equation (17)). As there is no straightforward Lyapunov function for the system, we monitor the number of unsatisfied clauses in real-time using a custom cost function $\kappa(\phi)$ defined as:

$$f_B = C_1 \wedge C_2 \wedge \dots \wedge C_{M-1} \wedge C_M \\ \rightarrow \kappa(\phi) = K_1 + K_2 + \dots + K_{M-1} + K_M$$

where $K_m(\phi) = l_m^m l_2^m l_3^m$ with literals $l_m^m = 0.5(1 \pm \tanh(\beta \cos \phi_j^m)) \in [0, 1]$. The sign that weights the tanh term depends on whether the literal $l_m^m$ corresponds to a positive $x_j(-)$ or negated variable $\bar{x}_j(+)$. This way, $K_m(\phi) = 0$ if there is a variable assignment that satisfies the clause $C_m$. Consequently, $f_B$ is true if $\kappa(\phi) = 0$. The tanh function is used to map phases to Ising spin values and is equivalent to rounding phases to the nearest multiples of $\pi$. For a sufficiently high $\beta$ value, we then have $\kappa(\phi) = N_{unsat}$, i.e. $\kappa(\phi)$ counts the number of unsatisfied clauses. Note that with the proposed rounding procedure, a clause $C_m$ is true if and only if $K_m(\phi) < 0.5^3 = 0.125$ ($K_m(\phi) = 0.125$ when $\forall j \cos \phi_j^m = 0$). Thus, if $\kappa(\phi) < 0.125$, $f_B$ is true and we
<figure><img src="image_13.png" /><figcaption>Figure 7. Simulation for the satisfiable SATlib instance ‘cnf-100-01’ with N = 100 variables and M = 430 clauses. Here, we monitor the cost function and stop the simulation when the cost<0.125 or when the system reaches a fixed point. (a) Cost function comparison between the standard ONN and the Lagrange version. While the two systems produce a rapid cost decrease in about 1 oscillation cycle (with more than 35 satisfied clauses), Lagrange oscillators are then actively exploring the phase space when the standard ONN gets stuck into a local minimum. When the Lagrange ONN finds an optimal phase assignment at t = 150 oscillation cycles, we take a snapshot of all energy terms Zm. (b) Standard ONN energy snapshot. The black stars show the energy values at the snapshot time. Most of them are not settling to the target Zm = 0. (c) LagONN energy snapshot. When the simulation is stopped (cost < 0.125), Zm values are getting closer to the origin. By monitoring the cost in real-time, we do not need to wait for full convergence towards an optimal saddle point where all Zm = 0.</figcaption></figure>

use this value as a threshold to stop the LagONN search as shown in figure 7(a). for N = 100 variables and M = 430 clauses.

Figures 7(b) and (c) show the energy values Zm for each clause at the initialization (purple cross) and at the snapshot time t = 150 oscillations (black star) when LagONN finds an optimal solution. For the nominal ONN, the final fixed point does not satisfy Zm = 0 for many clauses. For LagONN, since we stop the simulation before convergence, most of the final Zm values are zero yet. In practice, one could have a standard Boolean circuit corresponding to the formula $f_B$ (with AND and OR gates) checking in real-time the number of satisfied clauses and sampling the phases when the cost reaches the target value $\kappa(\phi) < 0.125$.

## Appendix B. ODE solver for LagONN’s state equations

LagONN is challenging to simulate at a large scale due to the number of coupled differential equations scaling with the number of clauses. To best select a suitable ODE solver, we studied the stiffness of LagONN’s state equations via a custom ODE solver in Matlab, which has an adaptive time step. We chose Fehlberg’s method, which has a local error scaling as O(dt⁴), requiring three evaluations of $\nabla L_T(\phi)$ per time step, and is a good compromise between speed and accuracy. It consists of a predictor/corrector method that provides a local error estimate at each time step, which is then used to adapt the step size. The predictor phase values $\phi_p$ are first computed according to Heun’s method as:

$$ \phi_p[k+1] = \phi[k] + dt[k] \times (f_1 + f_2)/2 \quad (21) $$

where $f_1 = \nabla L_T(\phi[k])$ and $f_2 = \nabla L_T(\phi[k] + dt[k] \times f_1)$.

The corrector phase values are then calculated according to Simpson’s rule as:

$$ \phi[k+1] = \phi[k] + dt[k] \times (f_1 + f_2 + 4f_3)/6 \quad (22) $$

where $f_3 = \nabla L_T(\phi[k] + dt[k] \times (f_1 + f_2)/4)$. We estimate the local error as:

$$ e_r = \sqrt{\frac{(\phi_p - \phi) \times (\phi_p - \phi)^T}{N + M}} \quad (23) $$

where N and M are the number of variables and clauses. Given a target error $\epsilon$ in radians per second, the new time step is calculated as:

$$ dt[k+1] = 0.9 dt[k] \Gamma \quad (24) $$

with $\Gamma = \sqrt{dt[k]\epsilon/e_r}$. If $\Gamma < 1$, the solver reiterates with a smaller time step. Otherwise, it integrates the next point with a larger time step. We varied the tolerated error $\epsilon \le 140^{\circ}$ using 10 instances per size N = 20 and N = 50 from the SATlib library [40]. Each LagONN instance was run 100 times with random initialization. Figure 8(a) shows the solver mean step size that exponentially increases with the tolerated error $\epsilon$.
<figure><img src="image_14.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 8.** (a) Median solver step versus the tolerated phase error $\epsilon$ for 10 3-SAT instances with $(N, M) \in (20, 91)$, (50, 218). Error bars correspond to 1st and 3rd quartiles. (b) 3-SAT cost vs tolerated error. (c) LagONN runtime vs tolerated error. (d) CPU runtime vs tolerated error.</figcaption></figure>

We found that LagONN is robust to numerical errors as it still finds optimal solutions with a similar runtime up to $\epsilon = 100^{\circ}$ (figures 8(b) and (c)). LagONN's runtime even decreases with numerical error, down to $10 \times$ for $N = 50$. However, when $\epsilon > 100^{\circ}$, i.e. $< dt > \approx 0.2$, LagONN's runtime significantly increases for $N = 50$. Based on this study, we ran our custom solver with a fixed time step $dt = 0.15$ for all simulations in the paper. The Matlab code was executed on a Linux server using one CPU per run.

## Appendix C. Algorithms

### C.1. LagONN

LagONN's pseudo code for Max-3-SAT is shown next, where we first construct the Lagrange function $L$ based on the input 3-SAT instance, and express its analytical derivatives used later on by the ODE solver for integration. In practice, both gradients of $L$ in the $\phi$ and $\phi_{\lambda}$ directions are computed by dedicated functions called by the solver 'integrate'. We used the custom ODE solver previously described in appendix B. We also construct a cost function $\kappa$ for the instance, which is called at each iteration to stop the dynamics when the phase assignment is optimal, as described in appendix A.

### C.2. Simulated annealing

In this paper, we execute SASAT, the simulated annealing algorithm proposed in [38] for SAT and expressed as follows:

The temperature decay rate was adjusted heuristically and scaled with the number of variables as $decay\_rate = 0.2/N$ to increase the annealing time and the success probability for larger instances. We have set $MAX\_TEMP = 1$ and $MIN\_TEMP = 0.01$ for all the experiments by inspecting the cost trace for several sizes.
**Algorithm 1.** Simulate LagONN for satisfiable 3-SAT.

**Input:** 3-SAT instance, Number of trials MAX_TRIALS, Simulation time MAX_TIME, solver time step dt

**Output:** A phase assignment $\phi$

1. Create Lagrange function $L$ for the instance;

2. Create gradient functions $\nabla_{\phi} L$ and $\nabla_{\phi_{\lambda}} L$;

3. Create cost function for the instance $\kappa$;

4. $N_{\text{iteration}} \leftarrow \lceil \text{MAX}\_\text{TIME}/dt \rceil$;

5. **for** trial $\leftarrow 1$ **to** MAX_TRIALS **do**

6. &nbsp;&nbsp;&nbsp;&nbsp;($\phi, \phi_{\lambda}$) $\leftarrow$ random phase;

7. &nbsp;&nbsp;&nbsp;&nbsp;**for** $i \leftarrow 1$ **to** $N_{\text{iteration}}$ **do**

8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;($\phi, \phi_{\lambda}$) $\leftarrow$ integrate($\phi, \phi_{\lambda}, \nabla_{\phi} L, \nabla_{\phi_{\lambda}} L, dt$);

9. &nbsp;&nbsp;&nbsp;&nbsp;**if** $\kappa(\phi) < 0.125$ **then**

10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** $\phi$;

**Algorithm 2.** SASAT algorithm [38].

**Input:** A set of clauses with $N$ variables, MAX_TRIALS, MAX_TEMP, and MIN_TEMP

**Output:** A variable assignment $S$

1. $\text{trial} \leftarrow 1$;

2. $\text{decay\_rate} \leftarrow 0.2/N$;

3. **while** trial $\le \text{MAX\_TRIALS}$ **do**

4. &nbsp;&nbsp;&nbsp;&nbsp;S $\leftarrow$ a random variable assignment;

5. &nbsp;&nbsp;&nbsp;&nbsp;j $\leftarrow$ 0;

6. &nbsp;&nbsp;&nbsp;&nbsp;T $\leftarrow \text{MAX\_TEMP}$;

7. &nbsp;&nbsp;&nbsp;&nbsp;**while** T $\ge \text{MIN\_TEMP}$ **do**

8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** S satisfies the clauses **then**

9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**return** S;

10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;T $\leftarrow \text{MAX\_TEMP} \cdot \exp(-j \cdot \text{decay\_rate})$;

11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $i \leftarrow 1$ **to** N **do**

12. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the cost change $\delta$ in the number of unsatisfied clauses if $i$ is flipped;

13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Flip $i$ with probability $1/(1 + \exp(\delta/T))$;

14. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;S $\leftarrow$ the new assignment if flipped;

15. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**j** $\leftarrow j + 1$;

16. &nbsp;&nbsp;&nbsp;&nbsp;trial $\leftarrow$ trial + 1;

# Appendix D. LagONN saddle points and dynamics

## D.1. Optimal saddle point

Here we prove theorem 2 which generalizes theorem 1 to M clauses.

**Theorem 2.** Let $L_T(\phi, \phi_\lambda) = \sum_m^M u_\lambda^m \vec{Z}_m$ and the formula $f_B = C_1 \wedge C_2 \wedge \dots \wedge C_{M-1} \wedge C_M$ is satisfiable, then:

1. $L_T$ has a at least one saddle point $L_T(\phi^*, \phi_\lambda^*) = 0$ such that $L_T(\phi^*, \phi_\lambda) \le L_T(\phi^*, \phi_\lambda^*) \le L_T(\phi, \phi_\lambda^*)$.

2. Such saddle point satisfies the constraints $Z_m(\phi^*) = 0$ for all clauses.

**Proof.** 1. We motivate the search for a saddle point using the concept of duality [67]. Consider the following dual function defined as:

$$D_T(\phi_\lambda) = \min_{\phi} L_T(\phi, \phi_\lambda). \qquad (25)$$

For any vector of phases $\phi_\lambda$ it is possible to find an optimal assignment of phase $\phi^*$ such that constraints are satisfied, i.e. for all clauses $Z_m = 0$ which gives as optimal value $L_T(\phi^*, \phi_\lambda) = 0$. Hence, $D_T(\phi_\lambda) \le 0$. The dual problem consists of finding the best lower bound for the optimal value of the initial problem—that is satisfying the constraints $Z_m = 0$ for all clauses $C_m$. Hence, we are looking for $\phi_\lambda$ that maximizes $D_T(\phi_\lambda)$ as:

$$\max_{\phi_\lambda} D_T(\phi_\lambda) = \max_{\phi_\lambda} \min_{\phi} L_T(\phi, \phi_\lambda). \qquad (26)$$

For any $\phi$, we can find a vector of phases $\phi_\lambda$ such that their corresponding unitary vectors are orthogonal to their corresponding $\vec{Z}_m$, hence, $\max_{\phi_\lambda} \min_{\phi} L_T(\phi, \phi_\lambda) = 0$. Consider now the inverse situation where
we first maximize $L_T$ as:

$$P_T(\phi) = \max_{\phi_{\lambda}} L_T(\phi, \phi_{\lambda}). \qquad (27)$$

For any $\phi$, we can set the phases $\phi_{\lambda}$ such that their corresponding unitary vectors point in the same direction as their $\vec{Z}_m$. Hence, $P_T(\phi) \ge 0$. Seeking the best higher bound is expressed as:

$$\min_{\phi} P_T(\phi) = \min_{\phi} \max_{\phi_{\lambda}} L_T(\phi, \phi_{\lambda}). \qquad (28)$$

For any $\phi_{\lambda}$, the best higher bound $P_T(\phi)$ is obtained when all the constraints are satisfied, i.e. $Z_m = 0$ for all clauses $C_m$. Hence, $\min_{\phi} \max_{\phi_{\lambda}} L_T(\phi, \phi_{\lambda}) = 0$. In summary, we obtain:

$$L(\phi^*, \phi_{\lambda}^*) = \max_{\phi_{\lambda}} \min_{\phi} L_T(\phi, \phi_{\lambda}) = \min_{\phi} \max_{\phi_{\lambda}} L_T(\phi, \phi_{\lambda}) = 0. \qquad (29)$$

Since $L_T$ is continuous in both $\phi$ and $\phi_{\lambda}$, and $L(\phi^*, \phi_{\lambda}^*)$ is attained for the satisfiable assignment of phase $\phi^* \in \{0; \pi\}^N$, $L(\phi^*, \phi_{\lambda}^*)$ is a saddle point satisfying the inequality $L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi, \phi_{\lambda}^*)$.

2. By contradiction, suppose that $(\phi^*, \phi_{\lambda}^*)$ is a saddle point of $L_T$ satisfying $L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi, \phi_{\lambda}^*)$ and $L_T(\phi^*, \phi_{\lambda}^*) > L_T(\phi^*, \phi_{\lambda}^*)$, but there is some clause $m$ such that $Z_m \ne 0$. Let us consider the corresponding Lagrange oscillator $u_{\lambda}^m$: the term $u_{\lambda}^m \cdot \vec{Z}_m$ is maximum when $u_{\lambda}^m$ points in the same direction as $\vec{Z}_m$. If at the saddle point $u_{\lambda}^m$ is not already pointing in the same direction as $\vec{Z}_m$, we can find a new Lagrange phase $\phi_{\lambda}$ for this clause such that $L_T(\phi^*, \phi_{\lambda}^*) > L_T(\phi^*, \phi_{\lambda}^*)$, which violates the saddle condition $L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi^*, \phi_{\lambda}^*)$ for all $\phi_{\lambda}$. If at the saddle point, the vector $u_{\lambda}^m$ with phase $\phi_{\lambda}^*$ points in the same direction as its respective $\vec{Z}_m \ne 0$, we can find a new assignment of phase $\phi$ satisfying all constraints and $Z_m = 0$ such that $L_T(\phi^*, \phi_{\lambda}^*) > L_T(\phi^*, \phi_{\lambda}^*)$, violating the saddle condition $L_T(\phi^*, \phi_{\lambda}^*) \le L_T(\phi, \phi_{\lambda}^*)$ for any $\phi$. $\square$

## D.2. Proposed dynamics to find a saddle point

To find an optimal saddle point as described by theorem 2, we combine gradient descent and ascent along $\phi$ and $\phi_{\lambda}$:

$$\left\{ \begin{array}{l} \tau \dot{\phi}_x = -\nabla_{\phi_x} L_T \\ \tau_{\lambda} \dot{\phi}_{\lambda} = +\nabla_{\phi_{\lambda}} L_T \end{array} \right. . \qquad (30)$$

Under these dynamics, the Lagrange function is not a Lyapunov function for the system since with our proposed dynamics $L_T$ can increase with time (gradient ascent along $\phi_{\lambda}$). Focusing on a single clause $i$, its time derivative is indeed expressed as:

$$\frac{dL_i}{dt} = \frac{du_{\lambda}}{dt} \cdot \vec{Z}_i + \frac{d\vec{Z}_i}{dt} \cdot u_{\lambda} \qquad (31)$$

and $L_i$'s gradient descent causes $\vec{Z}_i$ to evolve in the opposite direction from $u_{\lambda}$ as expressed here for $\tau = 1$:

$$\begin{aligned} \frac{d\vec{Z}_i}{dt} \cdot u_{\lambda} &= \frac{\partial \vec{Z}_i}{\partial \phi_x} \cdot u_{\lambda} \frac{d\phi_x}{dt} + \frac{\partial \vec{Z}_i}{\partial \phi_y} \cdot u_{\lambda} \frac{d\phi_y}{dt} + \frac{\partial \vec{Z}_i}{\partial \phi_z} \cdot u_{\lambda} \frac{d\phi_z}{dt} \\ &= - \left( \frac{\partial \vec{Z}_i}{\partial \phi_x} \cdot u_{\lambda} \right)^2 - \left( \frac{\partial \vec{Z}_i}{\partial \phi_y} \cdot u_{\lambda} \right)^2 - \left( \frac{\partial \vec{Z}_i}{\partial \phi_z} \cdot u_{\lambda} \right)^2 \\ &\le 0 \end{aligned} \qquad (32)$$

whereas $L_i$'s gradient ascent tends to bring $u_{\lambda}$ towards $\vec{Z}_i$ as:

$$\begin{aligned} \frac{d\vec{u}_{\lambda}}{dt} \cdot \vec{Z}_i &= \frac{\partial \vec{u}_{\lambda}}{\partial \phi_{\lambda}} \frac{d\phi_{\lambda}}{dt} \cdot \vec{Z}_i \\ &= (\vec{Z}_i \cdot \vec{u}_{\lambda})^2 \\ &\ge 0. \end{aligned} \qquad (33)$$
<figure><img src="image_18.png" /><figcaption>Figure 9. LagONN stability study for unsatisfiable 3-SAT instances (50 variables and 218 clauses from SATLib). The insets show the phase evolution for the N SAT variables. (a) LagONN running in its nominal mode with all oscillators activated. (b) LagONN with second harmonic injection locking (SHIL) to enforce phase binarization for the SAT variables. The SHIL amplitude is slowly ramped up until stabilization, similar to an annealing schedule. (c) LagONN simulation with the Lagrange oscillators stopped at t = 150 oscillations. The right plot shows the Lagrange phases $\phi_{\lambda}$ remaining constant after the stop.</figcaption></figure>

## Appendix E. LagONN stability for unsatisfiable problems

The stability of LagONN's dynamics is not guaranteed due to the gradient ascent mechanism introduced by the Lagrange oscillators. If the constrained problem is not satisfiable, such as an unsatisfiable 3-SAT formula, the oscillators never settle, as shown in figure 9(a) for a 50-variable unsatisfiable instance from SATLib. One way of stabilizing the system is to stop the Lagrange oscillator evolution for all clauses $m$ as $d\phi_{\lambda}^m/dt = 0$, for instance by setting the synaptic amplitude to 0 corresponding to a Lagrange time constant $\tau_{\lambda} \to +\infty$. In that case, the network has a Lyapunov function $L$ which is the Lagrange function itself. $L$ decreases over time as

$$ \begin{aligned} \frac{d}{dt} L(\phi, \phi_{\lambda}) &= \sum_{j} \frac{\partial L}{\partial \phi_j} \frac{d\phi_j}{dt} + \sum_{m} \frac{\partial L}{\partial \phi_{\lambda}^m} \frac{d\phi_{\lambda}^m}{dt} \\ &= -\tau \sum_{j} \left( \frac{d\phi_j}{dt} \right)^2 \le 0 \end{aligned} \quad (34) $$

which induces global stability as $L$ is bounded from below. An example of dynamics when stopping the Lagrange evolution at $t = 150$ oscillations is shown in figure 9(c), where the system settles to the optimal Max-3-SAT solution (one remaining unsatisfied clause). Stopping the Lagrange oscillators after some time is different from starting the dynamics without Lagrange oscillators, as the Lagrange phases $\phi_{\lambda}$ evolve before getting frozen to their final value.

Another option to enforce stability is the use of second harmonic injection locking (SHIL) to phase-lock oscillators to binary values. Since this technique is used in practice to mitigate oscillator frequency variations and recover binary Ising spins [17, 77], slowly annealing the injection strength could further freeze phases in the end before read-out, as shown with the simulation in figure 9(b) where we added a potential function $V_{SHIL} = -K(t) \sum_{i}^{N} \cos(2\phi_i)$ to the Lagrange function. In practice, $K(t) \ge 0$ is the increasing amplitude of a harmonic signal injected at $2\omega_0$, with $\omega_0$ the mean oscillator frequency. This injection could help scale up the system to overcome variability and stability issues in physical implementations.
## Appendix F. Impact of Lagrange oscillator speed

Throughout the paper, we assumed that the Lagrange and the standard oscillators are equally fast, i.e. $\tau = \tau_{\lambda} = 1$ in equation (17). Here, we study how speeding up or slowing down the Lagrange oscillators affects the runtime of the whole network. In particular, we set $\tau = 1$ and vary $\tau_{\lambda}$ for the Lagrange oscillator and measure the resulting TTS (see equation (18)).

Here we use the first 20 instances from SATlib with 20 and 50 variables with 100 trials for each instance. For each $\tau_{\lambda}$, we set a maximum simulation time $t_{max}$. Next, we check whether each trial finds an optimal solution in the predetermined simulation time. Based on these results, we compute the median success probability $p_s$ for each $\tau_{\lambda}$. Combining the maximum simulation time and the success probability, we can quantify the median TTS for 20 and 50 variables for each $\tau_{\lambda}$ value. The results are summarized in figure 10.

<figure><img src="image_19.png" /><figcaption>Figure 10. Impact of different $\tau_{\lambda}$ on the time-to-solution of LagONN. Here a simple benchmark on the Max-3-SAT problem with 20 and 50 variables for 20 instances each and 100 trials per instance. a.) Maximum simulation time as a function of $\tau_{\lambda}$ for 20 and 50 variables. b.) Median success probability $p_s$ as a function of $\tau_{\lambda}$ for 20 and 50 variables. c.) Median time to solution as a function of $\tau_{\lambda}$ for 20 and 50 variables.</figcaption></figure>

For minimum TTS, there seems to be an optimal value for $\tau_{\lambda} \approx 1$, which means all oscillators should have the same speed. Although more rigorous analysis is needed, $\tau_{\lambda}/\tau = 1$ may be linked to the saddle geometry of the optimal point (theorem 2), where minimization and maximization are interchangeable and maximization does not need to be performed in an outer loop as a slower process.

## Appendix G. Discretization of LagONN dynamics

To lay the groundwork for a potential LagONN digital hardware implementation, here we study the effect of phase discretization on the dynamics. We discretize the phase interval from 0 to $2\pi$ with a fixed number of states $N_{states}$. For example for $N_{states} = 16$ the phases can take on values

$$ \phi = \frac{2\pi}{16}k, \qquad (35) $$

with k an integer from 0 to $N_{states} - 1 = 15$. We study $N_{states} = 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192$ and run LagONN simulations on the first 20 instances from SATLib for 20 variables and 91 clauses, and for 50 variables and 218 clauses. Each instance is given 100 trials with random initial phases. We used the same random initial phases per trial, except discretized to the corresponding $N_{states}$. The TTS metric, as defined in equation (18), is used to compare different discretization levels. The results are shown in figure 11. For the case of 50 variables with $N_{states} = 16$ or $N_{states} = 32$, the data points are excluded due to unstable behavior.
<figure><img src="image_20.png" /><figcaption>**Figure 11.** Performance scaling when discretizing the phase dynamics of LagONN to different number of states $N_{\text{states}}$. The median runtimes are shown with the 25% and 75% quantiles as error bars.</figcaption></figure>

As can be seen in figure 11, there is a strong increase in TTS below $N_{\text{states}} = 64$ for instances with 20 variables, after which the TTS flattens. For 50 variables, the same can be seen, except that the trend flattens after around $N_{\text{states}} = 128$. We can conclude from these results that for a digital implementation of LagONN, one should take care to discretize the phases to a sufficiently high number of states to obtain a stable system. For 20 variables, this will be around $N_{\text{states}} \ge 64$, while for 50 variables it is around $N_{\text{states}} \ge 128$. Although more data is needed, we hypothesize that a higher number of variables beyond 50 would also require a higher number of $N_{\text{states}}$.

## References

[1] Cook S A 1971 The complexity of theorem-proving procedures Proc. 3rd Annual ACM Symp. on Theory of Computing (STOC'71) (Association for Computing Machinery) pp 151–8

[2] Mohseni N, McMahon P L and Byrnes T 2022 Ising machines as hardware solvers of combinatorial optimization problems Nat. Rev. Phys. 4 363–79

[3] Hopfield J J 1982 Neural networks and physical systems with emergent collective computational abilities Proc. Natl Acad. Sci. 79 2554–8

[4] Tank D and Hopfield J 1986 Simple ‘neural’ optimization networks: an A/D converter, signal decision circuit and a linear programming circuit IEEE Trans. Circuits Syst. 33 533–41

[5] Hopfield J J 1984 Neurons with graded response have collective computational properties like those of two-state neurons Proc. Natl Acad. Sci. 81 3088–92

[6] Lucas A 2014 Ising formulations of many NP problems Front. Phys. 2 5

[7] Venturelli D, Mandrà S, Knysh S, O’Gorman B, Biswas R and Smelyanskiy V 2015 Quantum optimization of fully connected spin glasses Phys. Rev. X 5 031040

[8] Parizy M and Togawa N 2021 Analysis and acceleration of the quadratic knapsack problem on an Ising machine IEICE Trans. Fundam. Electron. Commun. Comput. Sci. E104 1526–35

[9] Cellini L, Macaluso A and Lombardi M 2024 QAL-BP: an augmented Lagrangian quantum approach for bin packing Sci. Rep. 14 5142

[10] Hestenes M R 1969 Multiplier and gradient methods J. Optim. Theory Appl. 4 303–20

[11] Powell M J D 1978 Algorithms for nonlinear constraints that use Lagrangian functions Math. Program. 14 224–48

[12] Vadlamani S K, Xiao T P and Yablonovitch E 2020 Physics successfully implements Lagrange multiplier optimization Proc. Natl Acad. Sci. 117 26639–50

[13] Nagamatu M and Yanaru T 1996 On the stability of lagrange programming neural networks for satisfiability problems of prepositional calculus Neurocomputing 13 119–33

[14] Kirkpatrick S, Gelatt C D and Vecchi M P 1983 Optimization by simulated annealing Science 220 671–80

[15] Todri-Sanial A, Delacour C, Abernot M and Sabo F 2024 Computing with oscillators from theoretical underpinnings to applications and demonstrators npj Unconven. Comput. 1 14

[16] Parihar A, Shukla N, Jerry M, Datta S and Raychowdhury A 2017 Vertex coloring of graphs via phase dynamics of coupled oscillatory networks Sci. Rep. 7 911

[17] Wang T, Wu L, Nobel P and Roychowdhury J 2021 Solving combinatorial optimisation problems using oscillator based ising machines Nat. Comput. 20 287–306
[18] Graber M and Hofmann K 2022 A versatile and adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle 2022 IEEE 35th Int. System-on-Chip Conf. (SOCQ) pp 1–6

[19] Graber M and Hofmann K 2023 A coupled oscillator network to solve combinatorial optimization problems with over 95% accuracy 2023 IEEE Int. Symp. on Circuits and Systems (ISCAS) (IEEE) pp 1–5

[20] Delacour C, Carapezzi S, Boschetto G, Abernot M, Gil T, Azemard N and Todri-Sanial A 2023 A mixed-signal oscillatory neural network for scalable analog computations in phase domain Neuromorph. Comput. Eng. 3 034004

[21] Bashar M K and Shukla N 2023 Designing Ising machines with higher order spin interactions and their application in solving combinatorial optimization Sci. Rep. 13 9558

[22] Mallick A, Bashar M K, Truesdell D S, Calhoun B H, Joshi S and Shukla N 2020 Using synchronized oscillators to compute the maximum independent set Nat. Commun. 11 4689

[23] Ahmed I, Chiu P-W, Moy W and Kim C H 2021 A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems IEEE J. Solid-State Circuits 56 2870–80

[24] Moy W, Ahmed I, Chiu P-W, Moy J, Sapatnekar S S and Kim C H 2022 A 1,968-node coupled ring oscillator circuit for combinatorial optimization problem solving Nat. Electron. 5 310–7

[25] Graber M and Hofmann K 2024 An integrated coupled oscillator network to solve optimization problems Commun. Eng. 3 116

[26] Grollier J, Querlioz D, Camsari K Y, Everschor-Sitte K, Fukami S and Stiles M D 2020 Neuromorphic spintronics Nat. Electron. 3 360–70

[27] Grimaldi A, Mazza L, Raimondo E, Tullo P, Rodrigues D, Camsari K Y, Crupi V, Carpentieri M, Puliafito V and Finocchio G 2023 Evaluating spintronics-compatible implementations of ising machines Phys. Rev. Appl. 20 024005

[28] Dutta S, Khanna A, Assoa A S, Paik H, Schlom D G, Toroczkai Z, Raychowdhury A and Datta S 2021 An Ising Hamiltonian solver based on coupled stochastic phase-transition nano-oscillators Nat. Electron. 4 502–12

[29] Kim H W, Jeon S, Kang H, Hong E, Kim N and Woo J 2023 Understanding rhythmic synchronization of oscillatory neural networks based on NbOx artificial neurons for edge detection IEEE Trans. Electron Devices 70 3031–6

[30] Maher O, Jiménez M, Delacour C, Harnack N, Núnez J, Avedillo M J, Linares-Barranco B, Todri-Sanial A, Indiveri G and Karg S 2024 A CMOS-compatible oscillation-based VO₂ Ising machine solver Nat. Commun. 15 3334

[31] Yun S-Y, Han J-K and Choi Y-K 2024 A nanoscale bistable resistor for an oscillatory neural network Nano Lett. 24 2751–7

[32] Sharma A, Burns M, Hahn A and Huang M 2023 Augmenting an electronic Ising machine to effectively solve boolean satisfiability Sci. Rep. 13 22858

[33] Su Y, Kim T T-H and Kim B 2023 A reconfigurable Ising machine for boolean satisfiability problems featuring many-body spin interactions 2023 IEEE Custom Integrated Circuits Conf. (CICC) pp 1–2

[34] Hizzani M, Heittmann A, Hutchinson G, Dobrynin D, Van Vaerenbergh T, Bhattacharya T, Renaudineau A, Strukov D and Strachan J P 2024 Memristor-based hardware and algorithms for higher-order hopfield optimization solver outperforming quadratic ising machines 2024 IEEE Int. Symp. on Circuits and Systems (ISCAS) pp 1–5

[35] Pedretti G et al 2025 Solving boolean satisfiability problems with resistive content addressable memories npj Unconven. Comput. 2 7

[36] Dikopoulos E, Hsu Y-T, Wormald L, Tang W, Zhang Z and Flynn M P 2025 25.1 a physics-inspired oscillator-based mixed-signal optimization engine for solving 50-variable 218-clause 3-SAT problems with 100% solvability and 31.7µs solution time 2025 IEEE Int. Solid-State Circuits Conf. (ISSCC) vol 68 pp 01–03

[37] Salim A Y, Selman B, Kautz H, Ignjatovic Z and Köse S 2025 Ski-SAT: a cmos-compatible hardware for solving SAT problems IEEE Trans. Circuits Syst. I 72 1–12

[38] D Johnson and M Trick ed 1996 Cliques, Coloring and Satisfiability (Dimacs Series in Discrete Mathematics and Theoretical Computer Science) vol 26, (American Mathematical Society)

[39] Selman B, Kautz H A and Cohen B 1994 Noise strategies for improving local search Proc. 12th National Conf. on Artificial Intelligence (Vol. 1) (AAAI’94) (American Association for Artificial Intelligence) pp 337–43

[40] Hoos H 2000 Satlib - benchmark problems (available at: www.cs.ubc.ca/hoos/SATLIB/benchm.html) (Accessed 13 October 2023)

[41] Camsari K Y, Faria R, Sutton B M and Datta S 2017 Stochastic p-bits for invertible logic Phys. Rev. X 7 031014

[42] Cai F et al 2020 Power-efficient combinatorial optimization using intrinsic noise in memristor Hopfield neural networks Nat. Electron. 3 409–18

[43] Chowdhury S et al 2023 A full-stack view of probabilistic computing with p-bits: devices, architectures and algorithms IEEE J. Explor. Solid-State Comput. Devices Circuits 9 1–11

[44] Singh N S et al 2024 Cmos plus stochastic nanomagnets enabling heterogeneous computers for probabilistic inference and learning Nat. Commun. 15 2685

[45] Aadit N A, Grimaldi A, Carpentieri M, Theogarajan L, Martinis J M, Finocchio G and Camsari K Y 2022 Massively parallel probabilistic computing with sparse ising machines Nat. Electron. 5 460–8

[46] Aarts E and Korst J 1989 Simulated Annealing and Boltzmann Machines: a Stochastic Approach to Combinatorial Optimization and Neural Computing (Wiley)

[47] Fahimi Z, Mahmoodi M R, Nili H, Polishchuk V and Strukov D B 2021 Combinatorial optimization by weight annealing in memristive Hopfield networks Sci. Rep. 11 16383

[48] Jiang M, Shan K, He C and Li. C 2023 Efficient combinatorial optimization by quantum-inspired parallel annealing in analogue memristor crossbar Nat. Commun. 14 5927

[49] Goto H, Tatsumura K and Dixon A R 2019 Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems Sci. Adv. 5 eaav2372

[50] Ansari Md S, Nehete H, Grimaldi A, Mazza L, Raimondo E, Puliafito V, Finocchio G and Kaushik B K 2024 A comparison of oscillatory ising machines and simulated bifurcation machines for solving maximum cut problems 2024 IEEE 24th Int. Conf. on Nanotechnology (NANO) pp 389–93

[51] Ercsey-Ravasz M and Toroczkai Z 2011 Optimization hardness as transient chaos in an analog approach to constraint satisfaction Nat. Phys. 7 966–70

[52] Molnár B and Ercsey-Ravasz M 2013 Asymmetric continuous-time neural networks without local traps for solving constraint satisfaction problems PLoS One 8 e73400

[53] Traversa F L and Di Ventra M 2017 Polynomial-time solution of prime factorization and NP-complete problems with digital memcomputing machines Chaos 27 023107

[54] Suzuki H, Imura J-I, Horio Y and Aihara K 2013 Chaotic Boltzmann machines Sci. Rep. 3 1610

[55] Lee K, Chowdhury S and Camsari K Y 2025 Noise-augmented chaotic ising machines for combinatorial optimization and sampling Commun. Phys. 8 35
[56] Molnár B, Molnár F, Varga M, Toroczkai Z and Ercsey-Ravasz M 2018 A continuous-time MaxSAT solver with high analog performance *Nat. Commun.* **9** 4864

[57] Chang M, Yin X, Toroczkai Z, Hu X and Raychowdhury A 2022 An analog clock-free compute fabric base on continuous-time dynamical system for solving combinatorial optimization problems *2022 IEEE Custom Integrated Circuits Conf. (CICC)* pp 1–2

[58] Izhikevich E M and Kuramoto Y 2006 Weakly coupled oscillators *Encyclopedia of Mathematical Physics* (Academic) pp 448–53

[59] Jackson T, Pagliarini S and Pileggi L 2018 An oscillatory neural network with programmable resistive synapses in 28 nm CMOS *2018 IEEE Int. Conf. on Rebooting Computing (ICRC)* pp 1–7

[60] Abernot M, Gil T, Jiménez M, Núñez J, Avelido M J and Linares-Barranco B 2021 Théophile Gonos, Tanguy Hardelin and Aida Todri-Sanial. Digital implementation of oscillatory neural network for image recognition applications *Front. Neurosci.* **15** 1095

[61] Wang T, Leon W and Roychowdhury J 2019 New computational results and hardware prototypes for oscillator-based ising machines *Proc. 56th Annual Design Automation Conf. 2019 (DAC'19)* (Association for Computing Machinery)

[62] Burer S, Monteiro R and Zhang Y 2001 Rank-two relaxation heuristics for max-cut and other binary quadratic programs *SIAM J. Optim.* **12** 67

[63] Erementchouk M, Shukla A and Mazumder P 2022 On computational capabilities of ising machines based on nonlinear oscillators *Physica D* **437** 133334

[64] Bybee C, Kleyko D, Nikonov D E, Khosrowshahi A, Olshausen B A and Sommer F T 2023 Efficient optimization with higher-order ising machines *Nat. Commun.* **14** 6033

[65] Fisher M L 2004 The Lagrangian relaxation method for solving integer programming problems *Manage. Sci.* **50** 1861–71

[66] Wah B W and Wu Z 1999 The theory of discrete lagrange multipliers for nonlinear discrete optimization *Principles and Practice of Constraint Programming - CP'99 (Lecture Notes in Computer Science)* vol 1713, ed G Goos, J Hartmanis, J Van Leeuwen and J Jaffar (Springer) pp 28–42

[67] Boyd S P and Vandenberghe L 2023 *Convex Optimization* 29th edn (Cambridge University Press)

[68] Zhang S and Constantinides A G 1992 Lagrange programming neural networks *IEEE Trans. Circuits Syst.* **II** 39 441–52

[69] Corti E, Jimenez J A, Niang K M, Robertson J, Moselund K E, Gotsmann B, Ionescu A M and Karg S 2021 Coupled VO₂ oscillators circuit as analog first layer filter in convolutional neural networks *Front. Neurosci.* **15** 19

[70] Mitchell D, Selman B and Levesque H 1992 Hard and easy distributions of sat problems *Proc. 10th National Conf. on Artificial Intelligence (AAAI'92)* (AAAI Press) pp 459–65

[71] Rønnow T F, Wang Z, Job J, Boixo S, Isakov S V, Wecker D, Martinis J M, Lidar D A and Troyer M 2014 Defining and detecting quantum speedup *Science* **345** 420–4

[72] Hoos H H and Stützle T 2000 Local search algorithms for SAT: an empirical evaluation *J. Autom. Reason.* **24** 421–81

[73] Selman B, Levesque H and Mitchell D 1992 A new method for solving hard satisfiability problems *Proc. 10th National Conf. on Artificial Intelligence (AAAI'92)* (AAAI Press) pp 440–6

[74] Sajeeb M M H, Aadit N A, Chowdhury S, Wu T, Smith C, Chinmay D, Raut A, Camsari K Y, Delacour C and Srimani T 2025 Scalable connectivity for Ising machines: dense to sparse (arXiv:2503.01177)

[75] Pelofske E 2025 Comparing three generations of d-wave quantum annealers for minor embedded combinatorial optimization problems *Quantum Sci. Technol.* **10** 02

[76] Delacour C 2025 Self-adaptive Ising machines for constrained optimization (arXiv:2501.04971)

[77] Bashar M K, Mallick A and Shukla N 2021 Experimental investigation of the dynamics of coupled oscillators as Ising machines *IEEE Access* **9** 148184–90