Corentin Delacour

Stefania Carapezzi

Gabriele Boschetto

Madeleine Abernot

Thierry Gil

Aida Todri-Sanial

A Mixed-Signal Oscillatory Neural Network for Scalable Analog Computations in Phase Domain

Keywords: Oscillatory Neural Network, Analog Computing, Mixed-signal design, NP-hard problems Table 1: State-of-the-art ONN architectures Analog LC Digital Analog LC Analog LC Analog relaxation Analog relaxation (PTNO) Ring Oscillator Analog differential Analog relaxation Coupling Transformers Resistors Resistors Resistors Capacitors Capacitors Resistors B2B inverters Current sources with DACs Capacitors

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction 1.Oscillatory Neural Networks

The synchronization of oscillators is a fascinating phenomenon studied for a long time. In the XVII century, Huygens noticed that two identical clocks attached to the same beam synchronize to an anti-phase state where the two pendulums have the same frequency and move in opposite directions [START_REF] Huygens | Oeuvres complètes de Christiaan Huygens[END_REF]. It is in the 1950s that scientists imagined computing with coupled oscillators. The main motivation at that time was to replace the bulky and slow vacuum tubes used for signal to compute digital functions [START_REF] John Von Neumann | Non-linear capacitance or inductance switching, amplifying and memory devices[END_REF]. At the same time in Japan, Goto developed a similar paradigm called the parametron, which consists of an LC resonant circuit oscillating at one-half the driving frequency. Using transformers as coupling elements, Goto was able to build multiple large-scale digital computers with up to 9600 parametrons in 1958 [START_REF] Goto | The parametron, a digital computing element which utilizes parametric oscillation[END_REF]. However, parametrons became obsolete in the 1960s due to the advent of transistors that were faster and more scalable. [START_REF] Takahashi | A brief history of the japanese computer industry before 1985[END_REF].

With the emergence of Artificial Intelligence (AI) and neural networks, researchers brought back the idea of phase computing to solve complex tasks like pattern recognition [START_REF] Jackson | An oscillatory neural network with programmable resistive synapses in 28 nm cmos[END_REF][START_REF] Corti | Coupled vo2 oscillators circuit as analog first layer filter in convolutional neural networks[END_REF][START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] and NP-hard combinatorial optimization problems (COPs) [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF][START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF][START_REF] Mallick | Using synchronized oscillators to compute the maximum independent set[END_REF][START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF] that are challenging for conventional digital electronics. Inspired by Hopfield Neural Networks (HNN) [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] and by the Kuramoto model [START_REF] Izhikevich | Weakly coupled oscillators[END_REF], Aoyagi [START_REF] Aoyagi | Network of neural oscillators for retrieving phase information[END_REF], Hoppensteadt and Izhikevich [START_REF] Hoppensteadt | Pattern recognition via synchronization in phase-locked loop neural networks[END_REF] have proposed a computing paradigm called Oscillatory Neural Networks (ONN) that models the phase dynamics of coupled oscillators. ONNs are particularly interesting for solving COPs as they are dynamical systems that converge naturally to fixed points corresponding to the minima of some energy landscape. In other words, there is no external control that makes the neuron state evolve. Instead, all the phases evolve in parallel and in continuous time, enabling fast and energy-efficient inferences in the analog domain [START_REF] Delacour | Energy-performance assessment of oscillatory neural networks based on vo 2 devices for future edge ai computing[END_REF].

State-of-the-art ONNs

Recently, a new interest in ONN has risen thanks to the emergence of novel oscillating devices that enable the fabrication of efficient ONNs [START_REF] Csaba | Coupled oscillators for computing: A review and perspective[END_REF]. Such as, spin-torque and spin Hall devices [START_REF] Torrejon | Neuromorphic computing with nanoscale spintronic oscillators[END_REF][START_REF] Grollier | Neuromorphic spintronics[END_REF], micro-electromechanical systems [START_REF] Maffezzoni | Analysis and design of boolean associative memories made of resonant oscillator arrays[END_REF][START_REF] Hoppensteadt | Synchronization of mems resonators and mechanical neurocomputing[END_REF], and transition metal oxide devices are all candidates for implementing ONNs using their oscillatory behavior and synchronization properties [START_REF] Corti | Coupled vo2 oscillators circuit as analog first layer filter in convolutional neural networks[END_REF][START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF][START_REF] Shukla | Synchronized charge oscillations in correlated electron systems[END_REF][START_REF] Carapezzi | Advanced design methods from materials and devices to circuits for brain-inspired oscillatory neural networks for edge computing[END_REF][START_REF] Jackson | Oscillatory neural networks based on tmo nano-oscillators and multi-level rram cells[END_REF]. Beyond-CMOS devices are promising as they generally allow a compact oscillator design using a single device that could be scaled down to the nanoscale. Nevertheless, CMOS-based ONNs benefit from the mature CMOS technology which enables rapid ONN development and facilitates its co-integration with conventional digital circuits [START_REF] Jackson | An oscillatory neural network with programmable resistive synapses in 28 nm cmos[END_REF][START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF][START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF][START_REF] Graber | A versatile adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle[END_REF]. In this work, we focus on ONNs that compute in phase domain, i.e. with neurons that oscillate at the same frequency. However, note that it is also possible to compute with various frequencies [START_REF] Dutta | Spoken vowel classification using synchronization of phase transition nano-oscillators[END_REF][START_REF] Nikonov | Convolution inference via synchronization of a coupled cmos oscillator array[END_REF]. Regardless of the technology, we identify three important criteria for designing a competitive ONN that computes in phase domain. It should have:

1. Homogeneous oscillating frequencies 2. Compact and linearly programmable signed synapses 3. A scalable architecture Even with the mature CMOS technology, achieving perfect matching between hundreds of oscillators is unfeasible for small-scale oscillators due to device-to-device variations. Hopefully, some techniques can overcome frequency mismatches such as calibration or Sub Harmonic Injection Locking (SHIL). SHIL consists in driving the oscillators with a harmonic signal that can lock to a Fourier harmonic of the oscillating signal [START_REF] Neogy | Analysis and design of sub-harmonically injection locked oscillators[END_REF]. In case of large frequency mismatches, the injection of a strong SHIL signal ensures phase locking among the oscillators [START_REF] Bhansali | Gen-adler: The generalized adler's equation for injection locking analysis in oscillators[END_REF][START_REF] Todri-Sanial | How frequency injection locking can train oscillatory neural networks to compute in phase[END_REF]. The second criterion promotes synapses that are compact, programmable with signed weights, and have a value proportional to their conceptual weight. Some architectures can lead to a non-linear mapping between the conceptual weights and their hardware implementation [START_REF] Delacour | Mapping hebbian learning rules to coupling resistances for oscillatory neural networks[END_REF], or even be unknown due to the high complexity of the dynamics. Finally, the ONN architecture must be scalable to compete with conventional computing and solve large-scale problems involving thousands or millions of synapses. For this reason, we believe that the ONN architecture should be modular, i.e. to support the interconnection of smaller sub-ONNs to build a larger system and avoid the implementation of a fully-connected network.

Table 1 presents the state-of-the-art ONN architectures and their features. We only consider ONN computing in the phase domain and based on electrical oscillators. For solving COPs, a general approach is to map the input graph to the ONN where vertices are oscillators, and edges are synapses. Some architectures such as [START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF][START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF] are dedicated to finding the maximum cut of a graph with weights of the same sign, as the synapses only implement negative weights. The main drawback is that both coupling capacitors and resistors are required to program negative and positive weights, respectively. Other architectures using differential LC oscillators enable signed weights using resistors only [START_REF] Wang | New computational results and hardware prototypes for oscillator-based ising machines[END_REF][START_REF] Chou | Analog coupled oscillator based weighted ising machine[END_REF] but are not scalable on chip due to the bulky LC tanks and resistors. Digital ONNs are promising as they are scalable and modular, as demonstrated by Ahmed et. al. with their 560 ring oscillators chip [START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF]. A recent promising fully-analog architecture for solving COPs has also been proposed by Graber et. al. [START_REF] Graber | A versatile adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle[END_REF] that consists of 400 oscillators coupled with nearest neighbors.

A mixed-signal approach

In this work, we introduce a new mixed-signal ONN architecture, named Saturated Kuramoto Oscillatory Neural Network (SKONN) that leverages both analog and digital domains to satisfy the three ONN design criteria. SKONN takes inspiration from the state-of-the-art analog ONN architectures for which the dynamics evolve naturally in continuous time and can easily be described by phase models like Kuramoto's or Izhikevich's [START_REF] Izhikevich | Weakly coupled oscillators[END_REF], thus facilitating the exploration of potential applications. SKONN's main novelty consists in setting the computation and propagation in the analog and digital domains, respectively. Fig. 1 illustrates SKONN with 4 fully-coupled neurons. Digital propagation has several advantages such as greater noise immunity, a higher fan-out, and smoother interfacing with other digital circuits. Moreover, the separation between the computation and propagation induces a natural implementation of feed-forward synapses that have never been implemented in literature, to the best of our knowledge.

The paper is organized as follows. First, SKONN architecture is described and its dynamics are derived, highlighting the link between the Kuramoto model and SKONN, and its ability to solve NP-hard Weighted Max-cut problems. Then, we introduce a SKONN-PCB proof-of-concept that solves Weighted Max-cut problems with 9 nodes. A 16-neuron 65nm-ASIC chip recently taped out is also presented, demonstrating SKONN feed-forward ability with a XOR example. Finally, we report on SKONN's performance scaling in solving the Weighted Max-cut problem, and benchmark with state-of-the-art solvers on G-SET instances.

Methods

SKONN architecture overview

A mixed-signal oscillating neuron

A SKONN neuron consists of a relaxation oscillator producing analog and digital oscillations with period T at its input and output nodes, respectively. Fig. 2a shows the block diagram of the oscillating neuron. It consists of a hysteresis circuit that commands a shaper block to charge and discharge a capacitor C L with constant current I bias . The voltage across the capacitor V in i is fed back to the hysteresis comparator that switches between V DD and 0 when V in i reaches the thresholds V H and V L , thus producing oscillations. V out i holds the phase state in the digital domain, whereas V in i is the analog evolution of the oscillation. Note that the input impedance of the oscillator is purely capacitive in the ideal case so that any charge sent to the input node causes an instantaneous phase shift. The analog waveform V in i supports the computation and is separated from the digital propagation V out i , enabling a feed-forward propagation of the phase information.

Fig. 2c shows an example of feed-forward propagation between two oscillators. The computation occurs in the analog domain at the input node of neuron i that gathers the output signals from neuron j. The oscillator output signal is a square digital-like signal that carries the oscillator state and evolves until the phase dynamics settle to a fixed point. Choosing a triangular waveform at the analog input leads to simple yet rich phase dynamics that are similar to the Kuramoto model, which is known to have interesting computational properties [START_REF] Izhikevich | Weakly coupled oscillators[END_REF]. Moreover, it skips the use of bulky

V 1 in (t) V 2 in (t) V 3 in (t) V 4 in (t) V 1 out (t) V 2 out (t) V 3 out (t) V 4 out (
(W ij ) C ij Sign Amplitude Programmation V j out V i in sign(W ij ) C 0 C 1 C 2 C 3 (b)
SKONN synapse LC tanks needed for producing sinusoidal oscillations. The neuron voltage dynamics are expressed in Appendix A for completeness, although this work rather focuses on phase dynamics that are more suitable to study phase-based ONNs.

V j out V i in I ij I ij V i in V i out IN OUT Hysteresis Shaper V i out V i in ANALOG DIGITAL (a) SKONN neuron C L C L V j out 𝛿V =- Q ij C eq V i in I ij 𝚫𝛟 𝐢𝐣 = 𝛑 𝟐 𝚫𝛟 𝐢𝐣 = 𝛑 (c) SKONN principle V i out sign(W ij ) = -𝟏 C ij V j out V i in I ij IN OUT Neuron j C L IN OUT Neuron i V i out Q ij -Q ij V H V L 𝛿t 𝛿t

Synaptic design and weight sign

A SKONN synapse S ij consists of a capacitor C ij that transmits current pulses, i.e. charges Q ij , from the output of oscillator j to the input of oscillator i. C ij can easily be programmed using a capacitor bank and a register, as shown in Fig. 2b. Instead of propagating the sensitive analog signal, SKONN only transmits the oscillator phase information in a robust manner. The digital output voltage V out j is applied to C ij that creates current spikes holding the phase information ϕ j . The synaptic spike train can be expressed as follows:

I ij = C ij dV out j dt - dV in i dt (1)
The synaptic capacitor can be thought of as a digital-to-analog phase converter. The synaptic weight consists of the capacitance value C ij that linearly modulates the charge sent to the oscillating input node i as Q ij = C ij V DD , thus inducing phase shifts in the oscillation i as shown in Fig. 2c. To implement a negative weight, the complementary of

V out j defined as V out j (t) = V out j (t -T /2
) is selected using a multiplexer and applied to C ij . Compared to resistors, synaptic capacitors have several advantages for upscaling the ONN:

1. ONN computation models are generally based on the weak coupling assumption [START_REF] Izhikevich | Weakly coupled oscillators[END_REF][START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF] and necessitate weak synaptic signals. This means the ONN needs either large coupling resistors or small capacitors, the latter being much more scalable in a chip. 2. For a limited neuron output strength, the only way of increasing the synaptic fan-out is to reduce the synaptic current, which again would lead to bulky resistors or smaller capacitors in the case of SKONN.

SKONN phase dynamics 2.2.1 2 coupled oscillators

SKONN computing mechanism is illustrated in Fig. 2c with the case of a neuron j feeding its phase to another neuron i in a feed-forward manner with a negative weight. This means that V out j is selected by the multiplexer and applied to C ij , thus creating current spikes +Q ij δ(t) and -Q ij δ(t -T /2) that are aligned with the rising and falling edges of V out j , respectively. Each injected charge ±Q ij induces a voltage jump δV = ±Q ij /C eq at the input node, with The second phase catches up the first one after few cycles such that ∆ϕ ij ≈ 0 • . (c) Simulated dynamics with a negative weight W ij < 0. The second phase is pushed such that ∆ϕ ij ≈ 180 • after few cycles.

sign(W ij ) C ij MUX IN OUT V j out (t) V i in (t) V i out (t) C L I ij (t) (a) (b) IN OUT C L (c) W 21 > 0 W 21 < 0 ΔV = Q 21 C L Neuron j Neuron i V j in (t)
C eq = C L + C ij . As V in
i is a triangular waveform, δV provokes a time shift δt = ±C eq δV /I bias , where ± indicates here the sign of V in i 's slope. Knowing the period of the triangular oscillation T = 2C eq ∆V /I bias where ∆V is V in i 's peak-to-peak amplitude, we can then express the phase shift related to a single current spike:

δϕ = 2π δt T = π ±Q ij C eq ∆V = π ±C ij C eq V DD ∆V ≈ π ±C ij C L V DD ∆V if C L >> C ij (2)
SKONN's unique feature consists of this simple relation (2) between coupling capacitors and phase shift, thus enabling well-controlled phase dynamics and a precise weight mapping to the coupling capacitor C ij . As we will see later, the quantity β 0 = |δϕ/Q ij | provides the neuron phase sensitivity with respect to the charge perturbation. It is linked to the Phase Perturbation Vector (PPV) of the oscillator which is key for deriving SKONN's phase dynamics [START_REF] Bhansali | Gen-adler: The generalized adler's equation for injection locking analysis in oscillators[END_REF]. SKONN's PPV is defined and derived in Appendix B.

Two coupled oscillators converge either in-or out-of-phase, depending on the synaptic sign. To show this property, we use SKONN's phase dynamics that are derived in Appendix C using the PPV formalism [START_REF] Bhansali | Gen-adler: The generalized adler's equation for injection locking analysis in oscillators[END_REF]. Under the weak coupling assumption (C ij << C L ), the phase dynamics of oscillator i can be expressed as follows:

d dt ϕ i = 2β 0 Q ij T square(ϕ i -ϕ j ) (3)
With the 2π-periodic function

square(θ) = -1, if 0 < θ < π +1, if π < θ < 2π (4)
The phase fixed points can be derived from ( 3) and are expressed in the next proposition. Proposition 1. If the injected charge Q ij ̸ = 0 then the two SKONN oscillators admit a unique stable fixed-point

∆ϕ * = (ϕ i -ϕ j ) * such that ∆ϕ * = 0, if Q ij > 0 π, if Q ij < 0 ( 5
)
The proof is shown in Appendix C. In other words, propagating a spike train defined as ( 1) induces an in-phase or out-of-phase locking, depending on the polarity of Q ij . Each current spike produces a local phase shift to the analog input oscillation, resulting in an average phase shift ∆ϕ = ±2β 0 Q ij after each cycle (3). Fig. 3b shows a transistor-level simulation of the positive weight case. I ij perturbs V in i until the oscillators converge in phase. Similarly, Fig. 3c shows the same configuration with a negative weight and the oscillators are out-of-phase.

Note that the phases measured from the rising edges of V out i and V out j are slightly shifted from the theoretical fixed points [START_REF] Jackson | An oscillatory neural network with programmable resistive synapses in 28 nm cmos[END_REF]. This is mainly due to the limited bandwidth of the hysteresis block which does not switch instantaneously when reaching its thresholds. This non-ideality can be compensated and is further discussed in Appendix E. Interestingly, this phase shift disappears with symmetric synapses as both oscillators are equally delayed (see Fig. 6c).

N coupled oscillators

The phase dynamics of N sinusoidal coupled oscillators are often expressed using the Kuramoto model [START_REF] Izhikevich | Weakly coupled oscillators[END_REF][START_REF] Hoppensteadt | Pattern recognition via synchronization in phase-locked loop neural networks[END_REF][START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF]:

d dt ϕ i = -ω 0 N j=1 K ij sin ϕ i -ϕ j (6)
Where ω 0 is the frequency in rad/s and K ij the coupling coefficients. Similarly, we derive SKONN's phase dynamics for N oscillators as follows:

d dt ϕ i = ω 0 V DD ∆V N j=1 C ij C L square ϕ i -ϕ j (7)
Where we replaced β 0 and Q ij from (3) by their expressions

β 0 = π/(∆V C L ) and Q ij = C ij V DD . The derivation is detailed in Appendix C.
V DD is the digital voltage swing, ∆V is the peak-to-peak triangular voltage amplitude at the input, C ij is the synaptic capacitance value, and C L is the neuron input capacitance.

SKONN's phase dynamics are very similar to the Kuramoto model ( 6) except for its sinusoidal function replaced by a saturated square function in this work. It induces a binarization behavior that is useful for solving some optimization problems as shown later. Proposition 2. Consider a neuron i of degree D, i.e. driven by D neurons j with weighted charges Q ij ∈ {-q, +q} q ̸ = 0.

1. If D is odd and dϕ i /dt = 0, then there is at least one input neuron j such that (ϕ iϕ j ) is a multiple of π.

2. If D is even, then there is at least one ϕ i and one set of input phase ϕ j such that dϕ i /dt = 0 and ∀j (ϕ iϕ j ) is not a multiple of π.

The proof is shown in Appendix C. Interestingly, odd-degree neurons will phase-lock in-or out-of-phase with at least one input phase. The odd-degree property will be advantageous for solving some optimization problems on graphs,

Kuramoto SKONN

Even-degree neurons:

D = 2
Even and odd-degree neurons: as shown in the section 3.4. On another hand, an even number of inputs rather leads to a relaxed scenario where the neuron can settle into an infinite number of phases and could prevent phase binarization.

D 2 = 2 D 3 = 3
For instance, Fig. 4a shows a graph implemented with three fully-connected neurons with D=2. In this simulation, we randomly initialized the phases, numerically solved [START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] and measured the final phases with respect to the reference ϕ 1 . It can be seen from the distribution of final phases for 1000 trials that SKONN can settle to arbitrary analog phases, whereas the Kuramoto-ONN always converges to a single phase fixed point ϕ * = 0°, 120°, 240° . When SKONN has both odd and even numbers of input synapses, we heuristically find that most of phases tend to binarize as illustrated with the 5-node graph in Fig. 4b, although some phases (such as ϕ 2 ) can still converge to fixed points with arbitrary phase values. This aspect will be further discussed when solving larger graphs in section 3.4.

Note that similar dynamics have already been explored in simulation by Wang et. al. in their work about Oscillatory Ising Machines (OIM) [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF]. The authors studied the case where the sinusoidal term sin(∆ϕ) from Kuramoto ( 6) is replaced by tanh (α sin ∆ϕ) with α = 10. As SKONN's square interaction can be thought as square(∆ϕ) ≈tanh (α sin ∆ϕ) for α >> 1, we expect SKONN to have good performances when solving NP-hard combinatorial optimization problems.

SKONN energy landscape

SKONN stability can be proved by applying the Convergence Theorem for Oscillatory Neural Networks derived by Hoppensteadt and Izhikevich [START_REF] Frank | Weakly Connected Oscillators[END_REF]. With an odd coupling function (4) and symmetric coupling Q ij = Q ji , the theorem ensures that the phase differences converge to a stable fixed point. The proof consists of finding a Lyapunov function for the dynamics [START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] that is bounded below and minimized through time. A candidate for the SKONN Lyapunov function is:

E = β 0 T N i,j Q ij triangle ϕ i -ϕ j (8)
With:

triangle(θ) = θ -π/2, if 0 ≤ θ ≤ π 3π/2 -θ, if π ≤ θ ≤ 2π (9)
Under the assumption that Q ij = Q ji , one can check that:

∂E ∂ϕ k = β 0 T - N j=1 Q kj square ϕ k -ϕ j + N i=1 Q ik square ϕ i -ϕ k = - dϕ k dt (10)
Thus, SKONN minimizes E over time:

dE dt = N k=1 ∂E ∂ϕ k dϕ k dt = - N k=1 dϕ k dt 2 ≤ 0 (11)
SKONN's and Kuramoto's energy landscapes are represented for the three fully-coupled oscillators case in Fig. 5. We observe that SKONN's energy landscape has two plateaus where the phases can settle and remain stable, which is consistent with the phase distribution from Fig. 4a and linked to the even number of synaptic inputs (Proposition 2.2). Whereas the Kuramoto energy landscape consists of two minima and a saddle point, as highlighted by the three simulated trajectories of Fig. 5b andc. Unfortunately, we cannot visualise the energy landscape for larger networks but a local analysis around a phase fixed point can reveal the landscape around it (discussed in the supplementary material).

There can be sharp hills in some directions and plateaus in others. For instance in the 5-node case of Fig. 4b, most of the fixed points are surrounded by hills in all directions except for the direction ϕ 2 which consists of a trench where E remains constant.

So far, we have seen that SKONN has simple yet rich phase dynamics with unique phase binarization properties resumed in Proposition 2. Moreover, SKONN's phase evolution can be interpreted as the minimization of a N-dimensional energy landscape E (8) and thus naturally performs the gradient descent of E through time. One can harness this concept for solving optimization problems that consist in finding minima of a cost function. In the next sections, we focus on the NP-hard Max-cut problem and present two approaches for solving Max-cut with SKONN.

SKONN for solving the NP-hard Max-cut problem

Given a graph with a set V of N vertices connected by weighted edges W ij = W ji , the Max-cut problem consists in cutting the graph in two complementary subsets of vertices V 1 and V 2 such that the sum of weights between V 1 and V 2 is maximum. The Max-cut problem can be formulated as follows [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]:

Max 1 2 i,j W ij (1 -S i S j ) subject to: (12) S i ∈ {-1, +1} ∀ i ∈ V
Solving the Max-cut problem is NP-hard and the best-known approximation algorithm is the Semidefinite Programming (SDP) algorithm found by Goemans and Williamson [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] and denoted GW throughout the paper. By relaxing the binary spins S i to unit vectors v i in R N , GW relaxes the NP-hard Max-cut problem to an SDP convex problem for which optimality can be found in polynomial time:

Max 1 2 i,j W ij (1 -v i .v j )
subject to: (13)

v i ∈ R N |v i | = 1 ∀ i ∈ V
To compute the cut, the vectors are finally assigned to binary spins by splitting in two the N-dimensional sphere with a random hyper plan. Repeating this final rounding step provides a cut at least 0.878 times the maximum cut in expectation. However, due to the high dimension of the problem relaxation (R N ), GW is costly for large instances [START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF][START_REF] Honjo | 100,000spin coherent ising machine[END_REF] and alternative approaches using physical systems such as Quantum Annealers [START_REF] Hamerly | Experimental investigation of performance differences between coherent ising machines and a quantum annealer[END_REF], coherent Ising machines [START_REF] Honjo | 100,000spin coherent ising machine[END_REF], memristors [START_REF] Cai | Power-efficient combinatorial optimization using intrinsic noise in memristor hopfield neural networks[END_REF] or coupled oscillators are being investigated [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF][START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF][START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF][START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF].

The Ising approach

One of the most studied formalisms applied to ONN is from Ising which was initially derived to study magnetism in materials [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF]. Given interaction coefficients J ij ∈ R between particles that can have two spins S i ∈ {-1; +1}, the particles relax to a state that minimizes the Ising Hamiltonian (we skip the external fields for simplicity):

H = - 1 2 N i,j J ij S i S j (14)
Thanks to Lucas' seminal work [START_REF] Lucas | Ising formulations of many np problems[END_REF], all Karp's 21 NP-complete problems can be mapped to the Ising formalism and the solutions can be approximated by any physical machine that minimizes the Ising Hamiltonian [START_REF] Aoyagi | Network of neural oscillators for retrieving phase information[END_REF]. If SKONN phases take binary values ϕ i = (1 -S i )π/2 ∈ {0, π}, its Lyapunov function [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF] becomes:

E = β 0 T N i,j Q ij triangle( π 2 (S j -S i )) = - πβ 0 2T N i,j Q ij S i S j ∝ H (15)
Each synaptic current spike can be thought of as a downward step (due to equation ( 11)) in the energy landscape [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF] which corresponds to the Ising Hamiltonian (15) if the final phases are binary. However, having binary phases is not guaranteed in general as we have seen with the Proposition 2.2. To force phase binarization, it is common practice to inject into the oscillators a SHIL periodic signal at twice the oscillating frequency [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF] and described in Appendix D for SKONN.

The Max-cut problem can easily be mapped to an Oscillatory Ising Machine (OIM) with spins corresponding to the graph vertices by setting J ij = -W ij with W ij the graph weights [START_REF] Lucas | Ising formulations of many np problems[END_REF]. Then, the OIM performs the following minimization which is equivalent to the Max-cut [START_REF] Izhikevich | Weakly coupled oscillators[END_REF]:

Min H = Max - 1 2 i,j W ij S i S j subject to: (16) S i ∈ {-1, +1} ∀ i ∈ V
The general strategy is to 1) map the graph to the OIM, 2) start the OIM while ramping up a 2-SHIL signal to binarize the phases, and 3) read the stable phase state [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF]. Forcing phase binarization is common practice as it maps the OIM Lyapunov function to the Ising Hamiltonian [START_REF] Hoppensteadt | Pattern recognition via synchronization in phase-locked loop neural networks[END_REF]. However, how to binarize is not straightforward. If the injected signal is too strong, it may "freeze" the phases to sub-optimal local minima [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF]. Whereas if the signal is too weak, it might increase the OIM computation time. In this work, we rather harness the free SKONN dynamics without SHIL to compute the Max-cut.

A Rank-2 relaxation approach

Erementchouk et. al. [START_REF] Erementchouk | On computational capabilities of ising machines based on nonlinear oscillators[END_REF] have recently shown that the free OIM relaxation can be harnessed to solve the Max-cut problem. Their recent results demonstrate that letting a Kuramoto-ONN settle to analog phase values is equivalent to solving a rank-2 relaxation problem for the NP-hard Max-cut problem. Such phase dynamics are used in the CirCut solver [START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF]. Similarly to GW, the CirCut algorithm relaxes spins S i to 2D unit vectors x i ∈ R 2 such that x i = cos(ϕ i ) sin(ϕ i ) that can take arbitrary values on the unit circle. The objective of the rank-2 relaxation is:

Max 1 2 i,j W ij (1 -x i .x j ) = Max 1 2 i,j W ij (1 -cos(ϕ i -ϕ j ))
subject to: (17)

x i ∈ R 2 |x i | = 1 ∀ i ∈ V
Then, a rounding procedure produces spins to compute the graph cut. Unfortunately, this rank-2 algorithm cannot guarantee a lower bound on the cut as it remains a non-convex optimization problem. Nevertheless, its accuracy is comparable to the GW algorithm in practical use [START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF]. In this paper, we only explore the relaxation approach where we let SKONN settle without forcing binarization. We will see that SKONN's phase binarization property (Proposition 2.1) is particularly useful in this case. For the Ising approach, we invite the reader to consult the excellent work from Wang and Roychowdhury [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF] as the reported dynamics are equivalent to SKONN's.

Results

3x3 SKONN PCB

We designed a 3x3 SKONN on PCB with fully-connected capability and 81 synapses (Fig. 6a) as a proof of concept for the SKONN architecture. Due to area constraints, we only implemented negative weights that we program by placing discrete capacitors C ij manually. Fig. 6b shows the oscillating neuron based on a Schmitt trigger (U1) with feedback resistor R 3 that charges/discharges a load capacitor C L , producing a triangular-like waveform with 720 mVpp amplitude. The neuron output is a square-like waveform oscillating between V DD = +0.9V and V SS = -0.9V. Using an FPGA, we set the initial phase state by delaying the oscillator's starting time (Q1). The FPGA measures the neurons' output voltages and allows phase post-processing with a maximum precision of ϵ = 360 • f 0 /f F P GA . In our experiments we set f 0 = 4kHz, f F P GA =50 MHz and ϵ ≈ 0.03 • . Fig. 6c shows an experiment of two oscillators weakly coupled by C 12 = C 21 = 1%C L whereas Fig. 6d is a strong coupling with C 12 = C 21 = 10%C L . In both experiments, the oscillators are out-of-phase but the strong coupling case leads to a frequency reduction of -34% as the voltage jumps ∆V = (V DD -V SS )C 12 /C L produced by each current spike are too large with respect to the analog amplitude. This phenomenon can induce frequency mismatches between groups of strongly coupled oscillators and groups of weakly coupled oscillators. Frequency mismatches still need to be investigated and here we empirically choose C ij < 5%C L to guarantee phase locking among oscillators.

Weighted Max-Cut on PCB

The test case consists of the Max-cut problem with 2-bit positive weights. We generate random instances of Erdos-Rényi graphs G(N, p) [START_REF] Erdös | On the evolution of random graphs[END_REF] with N=9 nodes and p is the probability to have an edge between a pair of vertices such that the total number of edges m = p N (N -1)/2. For each graph edge, the weight is randomly selected from the list [0 10 22 47]/47 that corresponds to discrete capacitors used experimentally. Fig. 7a shows an example of a dense random graph instance with p = 0.75. We map the graph edges to the synaptic matrix and run 100 trials with random phase initializations. For each trial, the nine phases are sampled every oscillation period during 1000 oscillation cycles. Fig. 7b shows the final phases ϕ i (t = 1000 T ) measured for each trial, the latter indicated as the amplitude in the polar plot. The right-hand side of the polar plot corresponds to positive spins, whereas the left-hand side corresponds to negative spins. It appears that some phases such as ϕ 2 , ϕ 4 and ϕ 6 are always assigned to the same spin polarity whereas most of the phases can end up in both half-circles, depending on the phase initialization. Hence, SKONN final states depend on the initialization and several trials ensure obtaining a good solution. Fig. 7c shows the histogram of solution and the settling time. SKONN finds the graph Max-cut with 75% probability in less than 100 oscillation cycles on average. Fig. 8 presents another Weighted Max-cut problem with G(9, 0.5). Other Max-cut examples are shown in the supplementary material file.

4x4 SKONN CMOS design

SKONN Integrated Circuit

To further assess SKONN performances, we taped out an ASIC chip using a 65nm technology. Fig. 9a shows the chip layout view and matching to minimize the variations among oscillators and enhance the robustness. The second main objective was a low power consumption (160 µW without IOs). Although there are 16x more synapses than neurons, the synapses represent only 20% of the power consumption and are promising for solving dense networks. The constraint on the oscillator matching led to a large oscillator area (5000 µm²) but could be overcome with some calibration technique. Although the total chip area is 3.9 mm², there is room for improvement as the core area is 1.1 mm², and 71% of the total area consists of routing lines from the core to the high number of IOs and test pads (90).

Similarly to the SKONN design on PCB, the oscillating neuron consists of a hysteresis regenerative comparator whose digital output drives the charge and discharge of the capacitor C L =500fF in the shaper block (Fig. 9d). The comparator switches whenever the input voltage V in i reaches one of the two voltage thresholds V H and V L that define the analog voltage amplitude as ∆V = V H -V L =120 mV and (V H + V L )/2 = V REF =600 mV. The synaptic block consists of a capacitor bank ranging from 0 to 37.5fF that linearly maps the weight amplitude from 0 to 15

(C 0 = C 1 /2 = C 2 /4 = C 3 /8=2.5fF).
The weight sign is selected by the multiplexer commanded by the sign bit B S . The synaptic matrix is programmed by sending serially the 5x256=1280 bits through the registers.

Table 3 presents the specifications of the state-of-the-art ONNs designed using CMOS technology. The connectivity scheme can be all-to-all for small-sized ONNs (N≤100) but is obviously reduced to nearest neighbor connections for larger ONNs such as in [START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF] and [START_REF] Graber | A versatile adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle[END_REF] that chose to connect 6 and 8 (King's graph) neighbors, respectively. SKONN's input/output separation allows any type of modular connectivity in a robust manner as the synaptic signal from 2-coupled oscillators cannot leak to any other oscillator. This contrasts with fully-analog architectures such as [START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF][START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF] that merge input/output nodes, resulting in undesired current paths between non-adjacent oscillators. It appears that digital ONNs such as [START_REF] Jackson | An oscillatory neural network with programmable resistive synapses in 28 nm cmos[END_REF] and [START_REF] Ahmed | A probabilistic compute fabric based on coupled ring oscillators for solving combinatorial optimization problems[END_REF] are very energy-efficient as they produce a single oscillation with only 300 fJ. In contrast, analog oscillators found in [START_REF] Graber | A versatile adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle[END_REF] and this work consume 2.3 pJ and 10 pJ per oscillation, respectively. However, we believe that SKONN's energy could be reduced by relaxing the constraints on the analog oscillator matching and using a calibration scheme instead, as proposed by Graber et. al. [START_REF] Graber | A versatile adjustable 400 node cmos oscillator based ising machine to investigate and optimize the internal computing principle[END_REF].

Feed-Forward Network with SKONN

Propagating the information in a feed-forward manner is useful in some applications that require driving neurons. For instance, when training an ONN with the Equilibrium Propagation method [START_REF] Scellier | Equilibrium propagation: Bridging the gap between energy-based models and backpropagation[END_REF][START_REF] Zoppo | Equilibrium propagation and (memristor-based) oscillatory neural networks[END_REF], one must nudge the output

V i in (t) V i out (t) (b) (d) (c) (a) 2265 μm 1705 μm BS B0 B1 B2 B3 DIN DOUT LOAD CLK C L ΔV V j out (t) V i in (t)
ϕ XOR = 33°ϕ X = 160°ϕ Y = 185°ϕ Z = 0°X Y ϕ Z XOR -1 +1 +1 +1 +1 +1 +1 -1 Z H1 H2 ϕ Z W XZ W YZ ϕ Z (c) (a) (b)
W XZ W Y Z ϕ X ϕ Y ϕ XOR +1 +1 6°341°33°+ 1 -1 340°185°185°- 1 +1 185°340°185°- 1 -1 160°185°33°o
scillators toward the desired value which is challenging to obtain with recurrent synapses. Instead, teaching oscillators could drive the output oscillators using feed-forward connections without being impacted during the learning phase.

To demonstrate the SKONN feed-forward ability, we run ASIC transient simulations of a simple 2-input XOR operation. Inspired by the Parametron built by Goto in the 1950s [START_REF] Goto | The parametron, a digital computing element which utilizes parametric oscillation[END_REF], we use a 3-input SKONN neuron as a majority gate

ϕ M = (ϕ X .ϕ Y ) + (ϕ X .ϕ Z ) + (ϕ Y .ϕ Z ) (18)
Where ϕ X , ϕ Y , ϕ Z ∈ {0; 180} are the input binary phases thought as Boolean variables; i.e. ϕ M is true when ϕ M =180°. Interestingly, SKONN's odd-degree property (Proposition 2.1) ensures that ϕ M is binary when its inputs are also binary.

The XOR(X,Y) circuit is implemented by writing the XOR boolean expression ϕ

XOR = (ϕ X .ϕ Y ) + (ϕ Y .ϕ X ).
Considering ϕ Z as the reference oscillator, AND and OR gates are obtained when feeding ϕ Z or ϕ Z to the majority gates, respectively. Fig. 10a shows the obtained network with two hidden neurons (H1,H2), one neuron acting as a bias unit (Z), and one output neuron. The two input neurons X and Y are recurrently connected to the reference oscillator Z to set the input phase ϕ X and ϕ Y . When the weight W XZ is -1 or +1, it sets ϕ X ≈ 0 • or ϕ X ≈ 180 • , respectively, and the same rule applies for ϕ Y and W Y Z . All the oscillators are arbitrarily turned on at the same time before the input phases ϕ X and ϕ Y settle to the desired inputs. Then, the network further relaxes to a stable phase state after a few oscillations and we read the output phase ϕ XOR . Fig. 10b andc show the simulation results in the case ϕ XOR (180 • , 180 • ) = 0 • . Table 4 summarizes the results for the 4 possible inputs W XZ and W Y Z . By assigning the bit 0 when 270°≤ ϕ i ≤90°a nd 1 otherwise, it can be seen that the proposed network computes XOR(X,Y) in a feed-forward manner.

3.4 SKONN scaling and benchmarking

Weighted Max-cut of random graphs

To assess how SKONN's computational performances scale, we run large-scale simulations of random Weighted Max-cut problems for N=8, 16, 32, 64, 128, 256, 512, and 1024 nodes. For each graph density d=0.25, 0.5, and 0.75, 10 random graphs G(N, d) are generated such that the total number of edges m = d N (N -1)/2. The graph edges are randomly weighted with positive values from 0 to 15 that correspond to the ASIC synaptic range. We use the ASIC parameters and solve SKONN's dynamics [START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] with MATLAB using the built-in ODE solver ode15s. For each graph instance, we run 10 trials with random phase initialization, for a total of 100 trials per graph size and density. As a ground truth, we consider the best solution Cut GW provided by the Goemans-Williamson algorithm, out of 100 random projections defining the cut [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] and computed with the CVX solver on MATLAB [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]. The distance between the SDP cut Cut SDP and Cut GW is represented in Fig. 11a. As Cut SDP ≥Max-cut, the ratio Cut GW /Cut SDP gives a lower bound on the chosen GW cut. For all the trials, the results are compared with Kuramoto dynamics with the same parameters and initializations.

Fig. 11b shows SKONN's and Kuramoto's phase distributions for each ONN size. Here again, it appears that SKONN phases tend to be clustered near 0°and 180°whereas Kuramoto phases seem more uniformly distributed. Fig. 11c present the obtained cuts when considering the first oscillator as the reference, and normalized by Cut GW . We first notice that the results are quite homogeneous with respect to the graph densities. Secondly, it appears that SKONN produces high-quality cuts as Cut/Cut GW ≈ 1 for all ONN sizes. In contrast, Kuramoto-ONNs have a lower accuracy for sizes between N=16 and N=256. Interestingly, the settling time (time to reach a steady phase state) seems to grow according to a logarithmic law with the ONN size. This result refines some previous scaling observations mentioning a quasi-constant settling time [START_REF] Delacour | Energy-performance assessment of oscillatory neural networks based on vo 2 devices for future edge ai computing[END_REF][START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF]. It also confirms the high ONN parallelism and ability to compute in a few tens of cycles, even for large graphs.

Similarly to the CirCut algorithm [START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF] (rank-2 relaxation approach), we also investigate the Kuramoto accuracy when changing the reference oscillator and name it the "Kuramoto-CirCut" scheme. Fig. 11d presents the case where the best Kuramoto cut is chosen out of N possible cuts, whereas the SKONN reference oscillator remains the first one. It can be seen that SKONN provides the same quality cut as GW and Kuramoto-CirCut. However, compared to GW and Kuramoto-CirCut, SKONN's cut is solely obtained by reading out the phases with respect to a single oscillator and does not need N different cut evaluations that linearly increase the time to solution.

G-set benchmark

The previous study concerned random graphs. Here, we benchmark SKONN for solving Max-cut using the G-set benchmark that includes various graph topologies [START_REF] Martí | G-set data for max-cut[END_REF]. Table 5 shows the cuts obtained for a single trial with SKONN and Kuramoto using the same random phase initialization, and considering the first oscillator as the phase reference. The cuts are compared against the state-of-the-art GW [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF] and the Kuramoto-CirCut scheme [START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF]. GW's cut is the best cut obtained out of 100 random projections and computed with the CVX solver on MATLAB [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]. For graphs with N > 3000, GW values are taken from another state-of-the-art SDP solver [START_REF] Choi | Solving sparse semidefinite programs using the dual scaling algorithm with an iterative solver[END_REF] due to memory constraints. The Kuramoto-CirCut values correspond to the best cut extracted from the Kuramoto dynamics, out of N possible reference oscillators.

With a single run, Kuramoto-CirCut and SKONN solvers produce, on average, better results than the GW algorithm which picks up the best spin configuration out of 100 random projections. The average SKONN cut value is 94.6% of the best-known cuts [START_REF] Martí | G-set data for max-cut[END_REF] and the highest among the four methods. The accuracy obtained by simulating SKONN motivates its real hardware implementation as the time-to-solution could be drastically reduced compared to a CPU. For instance, solving the smallest graph G11 requires 12 s of GW runtime on a laptop (i7 Intel core @1.6 GHz and 32 GB of RAM). In contrast, SKONN's settling time does not vary much with the ONN size and could enable a large-scale cut computation in less than 431 cycles on average. With oscillators running at 1 MHz, the runtime per trial would only be 431 µs which is 2.8 10 4 x faster than GW's execution. However, reaching excellent cuts such as 99.9% of the best-known cut requires more trials and annealing the ONN to escape local minima. For more results using SHIL and various annealing schemes, please refer to [START_REF] Wang | Solving combinatorial optimisation problems using oscillator based ising machines[END_REF].

Discussion

SKONN mixed-signal architecture facilitates scaling up the ONN size thanks to the separation between analog computation and digital propagation. While the neurons exchange information in the digital domain, SKONN architecture ensures that the sensitive analog computation remains within modules and can be shielded from the outside. However, it remains some design challenges that are currently being investigated:

1. Any delay added to the digital propagation is equivalent to a synaptic phase shift. This must be considered for a high-frequency operation when buffers' delays become non-negligible and can be an issue for advanced technologies.

2. Large synaptic values can induce large voltage jumps at the analog input that can decrease the frequency for recurrent synapses, as highlighted in Fig. 6d. Without compensation, this strong coupling scenario could lead to multiple oscillation frequencies and prevent global phase-locking.

Conclusion

This article introduced SKONN, a novel mixed-signal ONN architecture that enables large-scale analog computations in phase domain. We presented how the association of simple circuitry can enable robust analog phase computation; while propagating the information in the digital domain to facilitate the implementation of large networks. We first reported experiments on a 9-neuron SKONN on PCB that finds the maximum cut of weighted graphs with high accuracy. Then, we presented the design of a CMOS 16-neuron integrated circuit (IC), highlighting SKONN programmability, modularity, and benchmarking with state-of-the-art ONN ICs. Furthermore, our design choice led to interesting phase dynamics that are a saturated version of the Kuramoto model and have unique properties. Such as, we found that such dynamics tend to binarize phases and are very efficient to solve NP-hard problems based on binary variables like Max-Cut. It appeared that SKONN's accuracy on random graphs is as good as the state-of-the-art Goemans-Williamson (GW) and CirCut algorithms, and even higher when benchmarked with graphs with up to 7000 nodes from the G-set.

Our study revealed that SKONN's computation time grows logarithmically with the network size which is promising for solving large-scale problems. For example, SKONN neurons oscillating at 1 MHz would provide a graph cut four orders of magnitude faster than the GW algorithm run on a CPU.

= V DD f H V in i , V out i , V L , V H -V out i (22)
With C eq the equivalent capacitance:

C eq = C L + j C ij (23)
(23) indicates that the synaptic capacitances are added to the oscillator load and slow down the charge and discharge of the input node. Large synaptic capacitances could potentially induce heterogeneous frequencies within SKONN and still needs to be explored.

B SKONN Phase Perturbation Vector

The PPV is a T-periodic function ⃗ v(t) that quantifies the phase shift of an oscillator subject to a perturbation occurring at time t [START_REF] Bhansali | Gen-adler: The generalized adler's equation for injection locking analysis in oscillators[END_REF]. One way of computing ⃗ v(t) is to inject a pulsed perturbation to the oscillator at time t, measure the induced phase shift and normalize by the perturbation's strength [START_REF] Dutta | An ising hamiltonian solver based on coupled stochastic phase-transition nano-oscillators[END_REF]. In SKONN, the synaptic current I ij perturbs the triangular oscillation V in i and the scalar PPV v(t) can be derived by computing the phase shift dϕ when injecting current pulses I(t ′ ) = dQδ(t ′t) with t ∈ [0; T [. From Fig. 13, we distinguish three cases: 1. 0 < t < T /2: the perturbed oscillation is shifted toward the left by the same amount of time -dt.

2. T /2 < t < T : the perturbed oscillation is shifted toward the right by the same amount of time +dt.

D Sub Harmonic Injection Locking in SKONN

To binarize the phases, one can inject a signal at twice the oscillating frequency V SHIL (t) = A sin(4πω 0 t) to an oscillator's node for which its scalar PPV contains a second-order harmonic P 2 ̸ = 0 in its Fourier decomposition [START_REF] Neogy | Analysis and design of sub-harmonically injection locked oscillators[END_REF].

In SKONN, we cannot inject the 2-SHIL signal to the input oscillating node V in i (t) as the associated scalar PPV only contains odd harmonics (square [START_REF] Hoppensteadt | Synchronization of mems resonators and mechanical neurocomputing[END_REF]). In practice, we inject the 2-SHIL signal to a biasing node that allows SKONN binary phase locking. In this case, the SKONN Lyapunov function becomes

E = β 0 T N i,j Q ij triangle ϕ i -ϕ j + N i A i P 2 cos(2ϕ i ) (50)
When SHIL amplitudes A i are large enough, phases are binarized ϕ i = (1 -S i )π/2 ∈ {0, π} and the SKONN Lyapunov function corresponds to the Ising Hamiltonian H with an additional offset:

E = - πβ 0 2T N i,j Q ij S i S j + N i A i P 2 = H + constant (51)

E Impact of SKONN's limited bandwidth

When SKONN is implemented in hardware, we observe some phase deviation with respect to the theoretical phase fixed points, as shown in Fig. 14. The main reason is the hysteresis block that does not switch instantaneously when the synaptic current spikes induce voltage jumps above or below the hysteresis thresholds V H and V L . To better understand this phenomenon, we ran two transistor-level simulations of two coupled neurons in feed-forward mode with a strong weight W 21 =+15 and different frequencies (Fig. 14). When the oscillation frequency is low (300 kHz), the hysteresis switching delay is negligible and there is only a small phase deviation δϕ=4 • . However, when the oscillation frequency increases to 1.2 MHz, the limited bandwidth of the hysteresis circuit causes a switching delay and a larger phase deviation δϕ=13 • .

For a given phase precision required by the application, this error could be mitigated by increasing the bandwidth of the hysteresis circuit or slowing down the oscillators, and constitute a trade-off with the energy consumption [START_REF] Delacour | Energy-performance assessment of oscillatory neural networks based on vo 2 devices for future edge ai computing[END_REF]. Interestingly, we have observed in experiments that having recurrent synapses W ij = W ji compensate the hysteresis delay induced in both neurons and the theoretical phase fixed point is reached in that case (see Fig. 6c andd).

MUX

Figure 1 :

 1 Figure 1: Illustration of a SKONN architecture composed of 4 oscillators and 16 synapses. Neuronal input and output lines are laid out vertically and horizontally, respectively. For neuron i, the input V in i (t) and output V out i (t) are synchronized such that V out i (t) commands the generation of V in i (t). A synapse S ij consists of a capacitor C ij that converts the digital signal V out j (t) in current spikes sent to the input node i. The multiplexer sets the weight sign by

  sign

Figure 2 :

 2 Figure 2: a) The SKONN neuron is a relaxation oscillator composed of two blocks: the hysteresis circuit that holds the neuron state and drives the shaper stage. The latter produces an analog triangular waveform at the input whereas the hysteresis block outputs a digital waveform. b) The SKONN synapse consists of a capacitor bank setting the weight amplitude. A multiplexer selects V out j or V out j to set the weight sign. c) SKONN principle of computation illustrated with a negative weight. The multiplexer selects V out j that is fed into C ij , creating current spikes I ij aligned with the rising and falling edges of V out j . The injected charges ±Q ij to V in i induce voltage jumps ±δV that cause time shifts δt. After a few cycles, the two oscillators lock to ∆ϕ ij = π.

21 21Figure 3 :

 213 Figure 3: (a) Illustration of a neuron that drives a second neuron in a feed-forward manner. The multiplexer sets the weight sign and the capacitor C ij converts the digital signal V out j (t) into current spikes I ij (t) that induce phase shifts of the input V in i (t). b) Transistor-level transient simulation with a positive weight W ij > 0.The second phase catches up the first one after few cycles such that ∆ϕ ij ≈ 0 • . (c) Simulated dynamics with a negative weight W ij < 0. The second phase is pushed such that ∆ϕ ij ≈ 180 • after few cycles.

Figure 4 :

 4 Figure 4: a) Three oscillators coupled by negative weights. The right-hand-side plots show the Kuramoto and SKONN distributions of the final phases for 1000 random initializations (uniform distribution). The Kuramoto-ONN converges to two fixed points (ϕ * 2 , ϕ * 3 )=(120°, 240°) or (ϕ * 2 , ϕ * 3 )=(240°, 120°). Whereas SKONN phases can settle to various analog phases as each oscillators has an even number of input synapses D. b) Five coupled oscillators where D 1 , D 3 and D 2 , D 4 , D 5 are odd and even, respectively. The Kuramoto-ONN settles to two fixed points (ϕ * 2 , ϕ * 3 , ϕ * 4 , ϕ * 5 ) ≈(248°, 137°, 331°, 166°) or (ϕ * 2 , ϕ * 3 , ϕ * 4 , ϕ * 5 ) ≈(112°, 223°, 29°, 194°). SKONN tends to binarize phases except ϕ 2 that can rather converge to any phase value.

Figure 5 :

 5 Figure 5: a) Three oscillators coupled by negative weights. b) and c) Corresponding Kuramoto and SKONN energy landscapes, respectively. The arrows represent three examples of trajectories with various initializations, highlighting Kuramoto's minima and saddle point. In this example, SKONN's energy minima consist of plateaus that are linked to the even-degree property (Proposition 2.2).

21 Figure 6 :

 216 Figure 6: (a) SKONN on PCB with 9 oscillators and 81 synapses. We set the weight amplitude with the synaptic capacitance value C ij . The FPGA initializes the phases by delaying the oscillators' starting time by switching Q1 and measures the digital oscillations buffered by the output stage (U2 and Q2). (b) The oscillator consists of a OP-amp Schmitt trigger fedback by the resistor R 3 to produce self-oscillations. R 1 //R 2 sets the analog oscillation amplitude. (c) Two coupled oscillators with C 12 /C L = C 21 /C L = 1%. (d) Two coupled oscillators with C 12 /C L = C 21 /C L = 10%.

Figure 7 :

 7 Figure 7: (a) Random instance of G(9, 0.75) with 2-bit weighted edges. (b) Phases measured after 1000 oscillation cycles. The polar amplitude represents the trial number (100 trials). We assign -90 • < ϕ i < 90 • → S i = +1 and 90 • ≤ ϕ i ≤ 270 • → S i = -1. (c) Histogram of Max-cut solutions for 100 trials and measured settling time. (d) Example of phase dynamics in polar representation. The polar amplitude corresponds to the time expressed in cycles. (e) Same dynamics represented in Cartesian plot with the cut evolution.

Figure 8 :

 8 Figure 8: (a) Random instance of G(9, 0.5) with 2-bit weighted edges. (b) Phases measured after 1000 oscillation cycles. The polar amplitude represents the trial number (100 trials). We assign -90 • < ϕ i < 90 • → S i = +1 and 90 • ≤ ϕ i ≤ 270 • → S i = -1. (c) Histogram of Max-cut solutions for 100 trials and measured settling time. (d) Example of phase dynamics in polar representation. The polar amplitude corresponds to the time expressed in cycles. (e) Same dynamics represented in Cartesian plot with the cut evolution.

Figure 9 :

 9 Figure 9: a) Layout view of the taped out SKONN ASIC. There are 16 neurons and 256 synapses laid out as in Fig.1. b) SKONN synapse with 5 bits precision. The four parallel capacitors set the the weight amplitude from 0 to 15. The bit B S selects the weight sign from the multiplexer. c) Shaper circuit producing the input triangular oscillation. d) Regenerative comparator that holds the neuron state and commands the shaper stage. An hysteresis behavior is obtained by setting the feedback transistors Q 7 and Q 8 such that W Q7 /W Q6 = W Q8 /W Q9 > 1 for the same transistor length.

Figure 10 :

 10 Figure 10: a) XOR(X,Y) circuit using SKONN and feed-forward synapses. Neuron Z is the reference oscillator and can be thought of as the neurons' bias, similar to perceptrons. The weights W XZ and W Y Z are the inputs and set the initial phases ϕ X and ϕ Y . b) Example of dynamics when W XZ = W XZ = -1 which set ϕ X = ϕ Y = 180 • and corresponds to the boolean case where X=Y=1. c) Zoom on the steady-state solution.16

Figure 11

 11 Figure11: a) Ratio between Cut GW and Cut SDP that gives a lower bound on Cut GW as Cut SDP ≥Max-cut. b) Phase distribution at steady state for SKONN (left) and Kuramoto (right) for various ONN sizes. c) Cuts obtained by SKONN and Kuramoto when the first oscillator is the phase reference and for various graph densities. d) Cut obtained when SKONN's reference is the first oscillator, compared to the case where the Kuramoto reference is changed N times, similarly to the CirCut algorithm[START_REF] Burer | Rank-two relaxation heuristics for max-cut and other binary quadratic programs[END_REF].

  π

Figure 13 :

 13 Figure 13: SKONN Phase Perturbation Vector (PPV). The injection of a charge dQ induces a time shift ±dt, which in turn creates a phase shift dϕ = 2π dt/T . The phase shift sign changes when there is a change in V in i 's slope.

3 .

 3 t ∈ {0; T /2}: the time shift is undefined as V in i 's slope is undefined (not derivable).Injecting dQ to C L induces a voltage jump dV = dQ/C L . This results in a time delay -dt and +dt when C L charges and discharges, respectively, with |dt| = C L dV /I bias . The amount of phase shift can then be expressed as dϕ = 2π dt/T . The oscillation period T is expressed by T = 2C L ∆V /I bias , with ∆V = V H -V L . Finally, merging the equations leads to dϕ dQ = ± π C L ∆V = ±β 0[START_REF] Maffezzoni | Analysis and design of boolean associative memories made of resonant oscillator arrays[END_REF]

Figure 14 :

 14 Figure 14: a) Two coupled oscillator in feed-forward mode with +C 21 =7.5%C L which implements W 21 =+15 in the ASIC. b) Transistor-level simulation with C L =2pF to decrease the frequency to f osc =300 kHz. c) Transistor level simulation in nominal case where C L =500fF and f osc =1.2 MHz.

  Table 2 contains the chip specifications. For this first ASIC version, we have focused on transistor

Table 2 :

 2 SKONN ASIC specifications

	Technology Neurons Synapses	Synaptic precision	f osc	SHIL inputs	Area	Power
	65nm	16	256	5 bits	1-4 MHz	16	3.9 mm² Neurons: 2 % Synapses: 24 % Biasing: 3 % Routing, IOs: 71 %	160 µW Neurons: 75 % Synapses: 20 % Biasing: 5 %

Table 3 :

 3 Comparison between state-of-the-art ONN Integrated Circuits

		Jackson et. al. [5] Ahmed et. al. [11] Bashar et. al. [19] Graber et. al. [20] This work
	Technology	28nm	65nm	65nm	28nm	65nm
	Neurons	100	560	30	400	16
	Connectivity	all-to-all	hexagonal	all-to-all	King's graph	all-to-all
	Power	303 mW	23 mW	1.76 mW	182 mW	160 µW
	Frequency	1 GHz	118 MHz	45 kHz	200 MHz	1 MHz
	Energy/osc	0.30 pJ	0.35 pJ	1.3 nJ	2.3 pJ	10 pJ
	Chip area	3.24 mm²	1.44 mm²	-	2.2 mm²	3.9 mm²

Table 4 :

 4 Solving XOR(X,Y) with SKONN

Table 5 :

 5 Comparison Between CirCut, Kuramoto, and SKONN solvers for G-SET Max-cut instances. The bold values indicate the best cut obtained among the four approaches. The number in parenthesis is the cut normalized by the best-known cut found in[START_REF] Martí | G-set data for max-cut[END_REF].

	Graph

Acknowledgments

We thank Laurent Deknyff from LIRMM, CNRS, for his circuit design expertise and for fabricating the PCB proofof-concept. We also thank Jeremie Salles and Fathi Ben Ali from LIRMM, CNRS, for their support during the IC design.

Declarations

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. Voltage (V)

Voltage dynamics of a single neuron obtained by solving numerically [START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF]. In this example, C L =500fF, I bias =200nA, V L =0.5V, V H =0.7V, V DD =1.2V, γ=100 and τ H =1ns.

A SKONN voltage and current dynamics A.1 Neuron voltage dynamics

By denoting I ij the input synaptic currents, the voltage dynamics of neuron i can be modeled as follows:

Where C L is the input capacitance, I bias the current that charges and discharges C L , V DD is the amplitude of V out i , V L and V H are the lower and upper thresholds of the hysteresis block, and τ H is the time constant linked to the output load of the hysteresis block. The term f H expresses the output switching with hysteresis behavior. As in [START_REF] Maffezzoni | Modeling and simulation of vanadium dioxide relaxation oscillators[END_REF], one can model the hysteresis behavior using a tanh function with slope γ:

=V DD and C L discharges. The two switching occur when V in i =V L and V in i =V H , respectively. Fig. 12 shows an example of numerical solution for the equations [START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF].

A.2 Synaptic currents

In SKONN, the digital output voltage V out j goes through the synaptic capacitance C ij that creates current spikes holding the phase information ϕ j . The synaptic spike train can be expressed as follows:

A.3 Voltage dynamics of coupled neurons

By injecting the synaptic current expression in equation [START_REF] Khairul Bashar | Experimental investigation of the dynamics of coupled oscillators as ising machines[END_REF], we obtain:

Which is the phase shift caused by the injection of 1 coulomb. The phase shift sign depends on whether the charge is injected during the charge or discharge of the triangular waveform. Considering the three previous cases and changing the time variable t to phase θ, we express SKONN's PPV as follows:

Where:

C SKONN phase dynamics

C.1 Two coupled oscillators

The phase dynamics of a single oscillator of frequency ω 0 = 2π/T are:

When the oscillator receives a pre-synaptic signal ⃗ b(t), it undergoes a time shift α(t) associated with the perturbation ⃗ b(t). If the variation of the oscillating amplitude remains small [START_REF] Bhansali | Gen-adler: The generalized adler's equation for injection locking analysis in oscillators[END_REF], α(t) dynamics can be expressed as follows:

Where ⃗ v(t) is the T-periodic Phase Perturbation Vector (PPV) associated with the oscillator; and describes the phase sensitivity of the oscillator under injections at different nodes. In our case, we consider scalars b(t) and v(t) as the pre-synaptic signal b(t) is injected to a unique input node. As b(t) also oscillates at frequency ω 0 and with phase ϕ b (t) = ω 0 t, we introduce ∆ϕ(t) = ϕ(t)ϕ b (t) = ω 0 α(t) that expresses the phase difference between post and presynaptic signals. The latter can be considered as the reference as it is driving the oscillator. To simplify equations, we define the 2π-periodic PPV and perturbation as v 2π (ω 0 t) = v(t) and b 2π (ω 0 t) = b(t), respectively. The phase dynamics become:

We assume that under weak coupling, the phase difference ∆ϕ evolves slowly compared to the presynaptic phase ϕ b and it is common practice to average out ∆ϕ over one period:

We saw previously that a pre-synaptic signal consists of current pulses that are aligned with the rising and falling edges of the digital pre-synaptic voltage. Then, we consider the case where b 2π (θ) consists of a train of Dirac pulses:

With:

) Under this assumption, (30) becomes:

In SKONN, the analog input oscillation is a symetric triangular waveform that has a simple 2π-periodic PPV expressed as follows:

And β 0 is the phase shift induced by the injection of 1 coulomb to the oscillating node. Finally, we express the phase dynamics of the driven oscillator:

We notice that the average of the phase dynamics are very similar to the Kuramoto model except that its sinusoidal interaction term is replaced by a saturated square function in our case. Proposition 1. If the injected charge Q ̸ = 0 then the two SKONN oscillators admit a unique stable fixed-point ∆ϕ * = (ϕ iϕ j ) * such that:

Proof. The proof consists in finding a Lyapunov function for the dynamics [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]. Consider the bounded continuous Lyapunov function:

With:

We have :

E is minimized through time as follows:

1. If Q > 0, the minima of E are ∆ϕ * = 0 [2π] and correspond to the phase fixed points of the dynamics (36).

2. If Q < 0, the minima of E are ∆ϕ * = π [2π] and correspond to the phase fixed points of the dynamics [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF].

In other words, propagating a spike train that consists of positive and negative current spikes spaced in time by T /2 induce an in-phase or out-of-phase locking, depending on the polarity of Q. The latter can be set by choosing one of the two complementary digital post synaptic voltages as shown in Fig 3.

C.2 N coupled oscillators

When an oscillator i is perturbed by N other oscillators with same pulsation ω 0 , (28) can be generalized in the scalar case:

Similarly to the two-oscillators case, averaging out the previous equation along the fast variable ϕ j leads to:

We use the spike train expression [START_REF] Neogy | Analysis and design of sub-harmonically injection locked oscillators[END_REF] to obtain:

As we inject pre-synaptic signals to the same node, we have v 2π ij = v 2π and we use the SKONN oscillator PPV v 2π (34) to finally get:

Note that we omitted the term ω 0 in the right-hand side of [START_REF] Scellier | Equilibrium propagation: Bridging the gap between energy-based models and backpropagation[END_REF] as in practice we refer to the relative phase relationship between oscillators instead of the absolute values that linearly increase with ω 0 t.

Considering SKONN's hardware implementation, we saw that β 0 is equal to:

And:

Where ∆V is the peak-to-peak triangular amplitude at the input, C L is the neuron input capacitance, V DD is the digital voltage swing, and C ij is the synaptic capacitance value. SKONN's phase dynamics become:

SKONN has an interesting phase binarization property resumed in the following proposition:

Proposition 2. Consider a neuron i of degree D, i.e. driven by D neurons j with weighted charges Q ij ∈ {-q, +q} q ̸ = 0.

1. If D is odd and dϕ i /dt = 0, then there is at least one input neuron j such that (ϕ iϕ j ) is a multiple of π.

2. If D is even, then there is at least one ϕ i and one set of input phase ϕ j such that dϕ i /dt = 0 and ∀j (ϕ iϕ j ) is not a multiple of π.

Proof.

1. By assuming the opposite, i.e. that ∀ j (ϕ iϕ j ) ̸ = 0 [π], it leads to ∀ j square(ϕ iϕ j ) = ±1, using (4). Noticing that D = m + l with m ̸ = l, and writing SKONN's phase dynamics [START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] leads to:

Which is not possible as m ̸ = l and proves the proposition.

2. Consider the integers m and l such that there are m weights Q ij = +q and l weights Q ij = -q, with D = m + l = 2k. We can choose:

(a) ϕ j =0 for the l (resp. m) input neurons. (b) ϕ j = π for kl (resp. km) other input neurons. (c) ϕ j =0 for the remaining k input neurons.

From SKONN's phase dynamics [START_REF] Delacour | Oscillatory neural networks for edge ai computing[END_REF] we obtain:

Assuming that ϕ i ̸ = 0 [π], it follows from (4):

T 2β 0 dϕ i dt = -ql(±1)q(kl)(±1) + kq(±1) = 0