Experimental Investigation of the Dynamics of Coupled Oscillators as Ising Machines

Coupled electronic oscillator networks, under second harmonic injection, have recently been shown to behave as Ising machines capable of solving computationally hard combinatorial optimization problems. In this work, we experimentally investigate the dynamical properties of a reconfigurable network of up to 30 oscillators (≡spins) configured as an Ising machine. Specifically, we analyze the characteristics of the solutions to the Ising model produced by the oscillators and show that as the system evolves towards the ground state through the high-dimensional phase space, it gets trapped in local minima resulting in sub-optimal solutions. Moreover, the exact local minima where the system gets trapped also changes implying that the trajectory of evolution of the system also changes with each trial. Finally, we illustrate experimentally how an appropriately designed annealing scheme can help the coupled oscillators escape a local minimum and attain a lower energy state.


I. INTRODUCTION
The Ising model, originally developed for spin glass systems [1], has recently received renewed attention for accelerating computationally hard problems that are considered intractable to solve using conventional digital computers.
The Ising Hamiltonian can be expressed as H = − ∑ , where (±1) is the i th spin, is the interaction coefficient between spin i and spin j, and N is the total number of spins in the system (the Zeeman energy term is ignored here). The goal is to find an optimal spin assignment that minimizes . Many computational problems such as computing the Maximum-Cut (MaxCut), Minimum Vertex Cover, Maximum Independent Set of a graph can be elegantly mapped to the Ising model [2]- [5]. For instance, consider the MaxCut problem which entails computing a cut that divides the nodes of the graph ( , ) into two sets ( 1 , 2 ) such that the number of common edges among the two sets is as large as possible. This problem can be mapped to by designating each node of the graph as a spin. When a node belongs to 1 ( 2 ), then it has a value (spin) +1 (-1), respectively; further, considering only binary weights (relevant to unweighted graphs), = −1, if an edge exists between nodes i and j; else = 0. This ensures that if the two adjacent nodes (i.e., nodes connected by an edge) are in the opposite sets, then = −1; if they are in the same set, = 1. Thus, computing the MaxCut of the graph entails minimizing ( =) − ∑ , which a physical Ising machine would accomplish by minimizing its energy to achieve the ground state. The applications of the Ising model have motivated an active effort to realize a technologically compatible physical Ising machine. Examples of existing approaches to realize such systems include quantum annealing systems (D-Wave) [6]- [10] which aim to overcome the fundamental complexity of the problem and facilitate an exponential speed up. However, quantum computing platforms require cryogenic operating temperatures [11] which implies that their application in energy-constrained environments (example, mobile applications, sensor networks) is limited. Besides quantum approaches, classical approaches have either relied on GPUbased & compute-in-memory-based digital implementations [12]- [17] that use various algorithmic methods to minimize the Ising Hamiltonian, or dynamical system-based designs such as, memristors [18]- [22], optical parametric oscillatorbased Coherent Ising machines [23]- [26], and coupled electronic oscillators [27]- [32] (considered here). In the latter approach, the physical system naturally aims to evolve towards the global energy minima, and subsequently, minimize the Ising Hamiltonian. The underlying motivation behind such designs is that while such classical systems may not overcome the theoretical hardness of the problem, the system will still offer better performance and higher energy efficiency during runtime for a large number of problems [4]. Moreover, such implementations do not require cryogenic cooling. Furthermore, among oscillator-based analog approaches, while optoelectronic oscillator-based CIM have shown promising performance, scaling up the CIM capacity while miniaturizing the hardware is challenging owing to its set up (which includes ~km long fiber) [33]. In contrast, electronic oscillators can be compact, integrated, and lowpower [34]- [44]; achieving these characteristics simultaneously in other systems may be potentially challenging [11], [33]. While the verdict on the relative computational performance of the various methods is still not out, coupled oscillators have shown early promise [45]. While our prior work [46], as well as other works [27]- [32], have focused on demonstrating the functionality of the oscillators as an Ising machine, this work differentiates itself by evaluating the dynamical properties of coupled oscillators and the characteristics of the resulting solutions to the Ising model. Specifically, we analyze the system evolution and the role of local minima in the phase space on the quality of the Ising solutions. We then experimentally show how annealing can assist the system in escaping such local minima.

A. COMPUTATION MODEL
The theoretical framework that describes how coupled oscillators under second harmonic injection behave as Ising machines and minimize the Ising Hamiltonian was developed by Wang et al. [45] and is summarized in this section. If an oscillator is perturbed by (or coupled with) another oscillator, then the generalized Adler's equation describing its phase dynamics can be written as, where, ( ) is the oscillator's phase, is perturbation's phase, is oscillator frequency, is perturbation frequency. is a periodic function known as the perturbation projection vector (PPV) and measures the response of an oscillator to a perturbation. For a system of N-coupled oscillators with a natural frequency , the generalized Adler's equation for the i th oscillator becomes, Considering that the oscillators are sinusoidal in nature (considered here for simplicity), oscillating at their natural frequency, and under the influence of second harmonic injection, equation (2) evolves to, Where, is the coupling coefficient between the i th and the j th oscillator (= −1, when an edge is present and = 0 when the edge is absent). The second harmonic injection helps discretize the phase to 0 and π [47]. Furthermore, it has been shown that there is a global Lyapunov function associated with the above equation, given by, It can be observed that equation (4) attains a minimum at the discrete points (0, π), which is equivalent to the Ising Hamiltonian with an offset. Further, this mapping is valid for non-sinusoidal oscillators [45].

B. COUPLED OSCILLATOR HARDWARE
Our experiments are performed using an integrated circuit (IC) consisting of 30 coupled oscillators with programmable capacitive coupling among them such that any one oscillator can be coupled to any and all other oscillators in the network. The IC is fabricated using bulk 65nm CMOS technology. Details of the IC have been described in our prior work [46] and are summarized below for completeness; additionally, the die photo of the oscillator IC and the PCB used for performing the measurements is shown in Fig. 1(a),(b), respectively, for reference. The oscillators are implemented using a CMOSbased Schmitt-trigger inverter where the oscillations are stabilized using a negative RC feedback ( Fig. 1(c)). A current mirror at the header and the footer of each oscillator helps regulate the oscillation frequency as well as facilitates the injection of the second harmonic signal. To enable the all-toall programmable coupling scheme described above, 870 (=P(30,2)) coupling elements are required with each element consisting of a coupling capacitor (CC=32.5fF) in series with a pass transistor (required to turn the coupling element ON & OFF) (Fig. 1(d)). The coupling scheme enables us to analyze a graph with arbitrary connectivity. Serial-input, paralleloutput (SIPO) shift registers are used to program the oscillators and the coupling network while the oscillator outputs are read serially using a multiplexer. The network to be analyzed is mapped to the oscillator hardware by creating a topologically equivalent oscillator circuit using the following relationship: spin ≡ oscillator, and spin interaction ≡ coupling capacitor CC (only binary interaction considered here), as shown in the example in Fig. 2(a). Fig. 2(a) shows the input spin network and the corresponding oscillator network onto which the spin network is mapped. The corresponding waveform ( Fig. 2(b)) of the coupled oscillators under the influence of the externally applied second harmonic signal, and the resulting phase plot (Fig. 2(c)) clearly shows a phase bipartition corresponding to the two spin states. The resulting spin assignment that gives rise to the minimum value of (= −3, here) is (1, 3) ↑; (2, 4) ↓.

A. STATISTICAL NATURE OF SOLUTION
While the coupled oscillator system evolves through the highdimensional solution space towards the ground state energy (corresponding to the optimal solution), it is likely to encounter many local minima where the system can get trapped, and subsequently, give rise to a sub-optimal solution. This is illustrated with the example of an oscillator network with 20 nodes and 114 edges (corresponding to the interactions among the oscillators) as shown in Fig. 3(a).

Equivalent
Oscillator Network  Figure 3(b) shows the measured solutions ( ) over 100 separate trials. Figure 3(c) compares the measured attained by the system (and its frequency) with the entire combinatorial solution space i.e., corresponding to all the possible spin assignments (grey in Fig. 3(c)) for the problem; there are 524,288 possible spin assignments out of which only 7 correspond to the optimal solution. It is evident from Fig. 3(c) that the spin configurations measured using the system are sub-optimal with the best solution only equal to 92% of the optimal value. Moreover, it can also be observed that the system energy (proportional to ) at which the peak of the distribution of the measured solutions (over the 100 trials) occurs, coincides with the value where the maximum number of local minima states occur. This indicates that in the absence of any external annealing, the system -despite trying to minimize its energy-gets trapped in one of the many local minima of the phase space, consequently, giving rise to suboptimal spin assignments. Furthermore, we also compute the Hamming distance among the measured spin assignments (Fig. 3(d)) to explore if there is a correlation among the solutions (generated in each trial). The resulting Hamming distance exhibits a Gaussian distribution implying that the solutions are widely different from each other and that the system gets trapped randomly in any one of the many local minima. This also indicates that the trajectory of the system in the phase space likely changes from run-to-run.  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication.
Where, =minimum value (optimal) of the Ising Hamiltonian, = maximum value of (=no. of edges in the network). is experimentally measured (expressed as ) using the coupled oscillators for various networks with 10, 20, and 30 oscillators, respectively; 20 randomly generated networks are considered for each case (i.e., graph size) across a broad range of network connectivity (5 graphs are considered for every ƞ= 0.2, 0.4, 0.6, 0.8 where ƞ is the edge density defined as the ratio of the number of edges in the network to the number of edges in an all-to-all connected network of the same size). 10 independent trials are performed for each network. It can be observed that while the oscillators compute a spin assignment in all the measured instances, the deviation of the measured solution from the optimal value increases with increasing network size; we observe that without annealing, the mean value of the measured solution is 96.4%, 84.5%, 83.7% of the optimal value for V=10, 20 and 30 nodes, respectively. This can be attributed to the increasing complexity of the solution space including the number of local minima having energy very close to the global minima. We also compute the Hamming distance between the solutions generated in the 10 trials for each graph (Fig. 4(d)) and observe a normal distribution for the Hamming distance which further corroborates the fact that the exact local minima where the system gets trapped is random and the system evolves with a different trajectory during every run.

B. ANNEALING
To prevent the system from getting trapped into a local minimum (resulting in sub-optimal spin assignments), we investigate the effect of using an annealing scheme. Fig. 5 illustrates the three different annealing schemes evaluated in this work; two (oscillator/spin) networks (G1, G2) with 30 nodes and randomly generated interactions (couplings) are tested for ƞ= 0.2, 0.4, 0.6, 0.8. While the first two schemes (A & B) modulate the amplitude of the injected second harmonic signal, the third scheme (C) turns the coupling ON and OFF, periodically. The modulation associated with the annealing scheme is applied three times and the best solution is selected. The corresponding reduction (Δ) in the mean value of (averaged over the 5 runs) before and after annealing is shown in Fig. 5.
Further, Fig. 6 compares the measured value of (represented using the normalized deviation d) for all the tested graphs (all 5 separate runs are considered) before and after annealing. It is evident that annealing helps the system achieve a lower energy configuration albeit not the optimal solution. The average improvement in the mean is 3.7% over all the measured configurations while the improvement in the best-case solutions (i.e., lowest ) obtained before and after annealing is 8.6%. Further, with annealing, the best solution (obtained within the 5 trials) is within 90% of the optimal solution in 6 (of the 8) measured graphs. While the number of networks tested is statistically insufficient to make a rigorous quantitative assessment of the improvement in the solution or to determine which of the annealing schemes is ideal, it is evident that incorporating an annealing scheme helps improve the performance of the system in all the experimentally measured cases.

IV. DISCUSSION
In summary, we have experimentally investigated the dynamical properties of coupled oscillator-based Ising machines using a testbed of 30 reconfigurable coupled oscillators and analyzed the evolution of the system through the phase space. Our results provide insights into the functional performance characteristics of Ising machines that would eventually be consequential to their practical deployment. We reveal the impact of the local minima on degrading the quality of solutions produced by electronic oscillator-based Ising machines, and subsequently, illustrate how an appropriately designed annealing scheme can improve the solution. Consequently, our work inspires the exploration of new pathways to optimize the performance of oscillator Ising machines. For instance, engineering an appropriate annealing schedule for a given problem could possibly leverage machine learning techniques. Additionally, hybrid digital-analog approaches, where digital post-processing algorithms are designed to help improve the solution quality of the oscillator Ising machine with minimal time penalty, can also be explored.