Analog Coupled Oscillator Based Weighted Ising Machine

We report on an analog computing system with coupled non-linear oscillators which is capable of solving complex combinatorial optimization problems using the weighted Ising model. The circuit is composed of a fully-connected 4-node LC oscillator network with low-cost electronic components and compatible with traditional integrated circuit technologies. We present the theoretical modeling, experimental characterization, and statistical analysis our system, demonstrating single-run ground state accuracies of 98% on randomized MAX-CUT problem sets with binary weights and 84% with 5-bit weight resolutions. Solutions are obtained within 5 oscillator cycles, and the time-to-solution has been demonstrated to scale directly with oscillator frequency. We present scaling analysis which suggests that large coupled oscillator networks may be used to solve computationally intensive problems faster and more efficiently than conventional algorithms. The proof-of-concept system presented here provides the foundation for realizing such larger scale systems using existing hardware technologies and could pave the way towards an entirely novel computing paradigm.


Introduction
Solving certain classes of combinatorial optimization (CO) problems has proven to be notoriously difficult using standard von Neumann computing architectures.A canonical example is the traveling salesman problem, for which exact algorithms scale very poorly with problem size.Applications of CO problems span many disciplines, including business operations, scheduling, traffic routing, finance, big data, drug discovery, machine learning, and many other systems requiring the minimization of a complex energy landscape with multivariate inputs 1 .As the decades-long progress of digital CMOS technologies begins to plateau, there is a growing desire to find alternative computing methodologies which can address the challenges of these traditionally difficult problems.Novel quantum annealing machines 2 and digital CMOS annealing accelerators are garnering significant attention in the hopes of solving such problems faster than conventional algorithms performing heuristic optimizations. 2,3One strategy to solve CO problems lies in mapping them to the Ising Hamiltonian for spin glass systems and finding the ground state solution 4,5 .The Ising model dates back many decades but was re-popularized by D-Wave Systems in an attempt to exploit quantum mechanical phenomena to speed up computation.The Ising Hamiltonian can also be mapped to the comparable field of Hopfield neural networks, which have been previously explored to solve CO problems 1,[6][7][8] .
Recently, alternative classical methods to solve the Ising model have emerged using optoelectronic parametric oscillators [9][10][11] , memristor cross-bar arrays 8,12,13 , electronic oscillators 14 , and GPU based algorithms 15,16 .Additionally, an analysis of one such coupled oscillator system revealed the potential for a significant speedup over digital computing algorithms at large node sizes 17 .Even algorithmic approaches which emulate nonlinear dynamical systems but are run on conventional computing hardware have been shown to match and even surpass the performance of state-of-the-art algorithms, motivating the desire to build these systems in physical hardware 18,19 .However, scaling up the optoelectronic oscillator Ising machine remains challenging owing to its time-multiplexed architecture and highly complex and costly optoelectronic setup 9,10 .The all-electronic oscillator concept initially proposed by Wang and Roychowdhury introduces the tantalizing prospect of creating a similar system using readily available electronic components interconnected in a parallel fashion and is particularly well suited for chip-scale integration and scaling using present day technologies. 14 this paper, we build upon this initial work and demonstrate a 4 node, fully connected LC oscillator based analog circuit with standard electronic components which accurately maps to the Ising model.We provide a thorough analysis of system performance, describe in detail the circuit architecture including a new modality used to program variable interconnection strengths, and demonstrate the system capabilities in solving a variety of MAX-CUT problems with both binary and multi-bit weight values.To the best of our knowledge, this is the first demonstration of an all-electronic oscillator-based weighted Ising machine.The statistical analysis presented herein provides valuable insight into the viability of these systems as computing platforms when scaled to larger node counts.

Theory
Much of the theoretical framework for using coupled oscillators as Ising solvers has been described previously, but will be summarized here for convenience. 20The Ising Hamiltonian is given by the following equation: where  represents the number of nodes in a particular problem set,   represents the weight values interconnecting the nodes, and  = [  …   ] represents the solution space where   can take the value of either +1 (spin ↑) or -1 (spin ↓).One common benchmark optimization problem, known as the MAX-CUT problem, is defined by taking an undirected graph and finding a bisection of that graph which maximizes the cut set.This problem can be mapped directly to the Ising Hamiltonian above, and the solution to the problem is represented by the state s which minimizes H.A particular problem of interest is typically defined by a graph G(V,E), where V represents the number of vertices and E represents the number of edges.An example of a 4-node system is shown in Fig. 1(a).If  = 1 for all connections, and we neglect the Zeeman term (ℎ  = 0), the minimum solution is defined by 6 degenerate 2x2 solution sets, as shown in Fig. 1(b).It has been shown previously that these graphs can be represented by a network of coupled nonlinear oscillators whose phase dynamics are described by the Kuramoto model, and that this model maps directly to the Ising Hamiltonian if the phases of these oscillators take values of either 0° or 180°1 4,21 .One mechanism used to binarize the phases is to introduce an injection-locking signal at twice the natural frequency of the oscillators 14 .This 'super-harmonic' injection locking signal, and its application in the context of LC-oscillator systems is mathematically similar to the case of the degenerate optical parametric oscillator used in previous optical Ising machines 9 , where optical pump pulses at twice the optical oscillation frequency are used to create binary phase values.In the LC-oscillator system, the phase evolution of each oscillator in the network (  (t)) can be represented by the following differential equation: Where N is the total number of oscillators (or vertices) and A represents the amplitude of the injectionlocking signal applied to each oscillator.The time t in the equation is defined in dimensionless units in relation to the oscillation period.The addition of Gaussian phase fluctuations, representing noise in the system, converts equation Eq. ( 2) into a network of stochastic differential equations (SDEs), whose solutions can be approximated iteratively using the Euler-Maruyama method.Annealing is accomplished by gradually increasing the injection locking term A according to the formula: with τ set to 5 oscillation cycles.One solution simulated for the example problem described above is shown in Fig. 1(c), demonstrating that starting from randomized initial conditions, the oscillator phases settle to one of the 6 correct solutions.The system settles to the ground state within 6 oscillation cycles.When this simulation was run 1000 times, the system always settles to a correct solution state with equal probability Fig. 1(d).

Results a. Circuit Implementation
The oscillator circuit employs a differential injection-locked frequency divider topology, as shown in Fig. 2(a) 22 .This architecture is selected to offer the option of synchronizing the oscillator with an incident superharmonic injection locking signal.Transistors M1 and M2 form a cross-coupled pair, which serves as the negative resistance component, necessary for unity loop gain.The coupling signal from the other oscillators is applied differentially through transistors M3 and M4.Current source I2 provide the bias current for the coupling signal.The output voltage of the oscillator is tapped at nodes VoL and VoR, directly out from the oscillating LC tank.The LC tank circuit has a resonant frequency of 50 kHz, and is composed of L1,C1, and L2,C2.Current source I1 provide a biasing current for the oscillator circuit.The source I1 can also be used for injection locking to help polarize the phases to 0º and 180º.
A differential analog multiply and accumulate circuit is used apply the necessary signals into input nodes ViL and ViR, as shown in Where  =   * 1024/20kΩ, and  is the mapping scalar.To prevent high currents at the summing amplifier at high coupling coefficient values, we set   = 760 .The mapping from the coupling coefficients,   , to the resistance values is scaled based on the  term to enable maximum dynamic range.A value of  = 2.5 is used in this paper.

b. Experimental Setup
A schematic of the experimental setup for the 4 node system is shown in Fig. 3  A representative phase-time trace of a single run is shown in Fig. 3(b).The digital potentiometers are programmed sequentially in time, due to the serial nature of the digital scheme.The time period of the weight programming is dependent on the loop speed of the Arduino board and the complexity of the code which has a maximum loop frequency of 117 kHz. 23After the last weight (J34) is programmed, the oscillators arrive at the new solution state.Phases of the oscillators are each measured relative to oscillator 1, whose phase is given a default binary value of "0".The binary phase outputs are then determined with a simple threshold, where phases differences > 90° with respect to oscillator 1 are assigned "1", and phases < 90º are assigned "0".We note that, unlike the simulation, no annealing of the coupling strength, is required to arrive at the solution.This is likely due to the low node count of the system.With more nodes, simulations of this oscillator network using coupled differential equations have shown that ramping up either the coupling strength or an external super-harmonic injection locking signal will be required in order for the system to preferentially settle to the correct solution.The super-harmonic injection locking signal serves a dual purpose in that it replaces the rounding feature described above, forcing the phases of each oscillator to their respective binary values (either 0º for "0" or 180º for "1").As pointed out previously, this signal is required to extract sensible answers using large node count systems. 14high sampling rate phase measurement of the phase transition of the oscillators is shown in Fig. 4. The phases were obtained via performing a moving window Fourier transform on the time-domain amplitude signal.The red dashed weight change line is when the last potentiometer (J34) has been digitally switched.
From this measurement, we observe a 100 µs solution time, which translates to 5 cycles of the 50 kHz oscillators, as shown in Fig. 4(a).To understand the relationship of the oscillator frequency to the timeto-solution, Fig. 4(b) shows the same experiment but with oscillator frequencies reduced by one order of magnitude, to 5 kHz, by utilizing 10 uF capacitors.We experimentally observe the time-to-solution also increases by approximately one order of magnitude, which demonstrates that the time-to-solution scales directly with the oscillator frequencies.The 90º threshold algorithm translates the final oscillator phase states to a binary solution of [0 1 1 0], if we use the order of the oscillator number.This solution is one of   10 V.With the appropriate binary weights shown in Fig. 5(a)-(d), each oscillator can be made to be 180° out of phase with all other oscillators.With this graph structure of only one node connected to the remaining three, it intuitively follows that the energy minimization is achieved when the central node is out of phase with the others.Under the fully connected system in Fig. 5(e), there exist 6 degenerate solutions, in which any 2 oscillators are 180º out of phase with any other 2 oscillators.Here we only show one of the six solutions.

c. Statistical Analysis
In order to validate the accuracy and reliability of the system, more than 2000 automated computer controlled experiments with random weights were performed.Fig. 6(a) and 6(b) shows the probability distribution of the arrival of each solution energy with a 1 bit and 5 bit weight values, respectively.The probability of achieving the ground state for 1 bit and 5 bit weights are 98% and 84%, respectively.The ground state probability was determined by 10 trials of each of the 2000 randomly generated problems.
Clearly, the probability of achieving the ground state decreases with increasing weight bit resolutions.We observe that the solutions for all of the experiments follow a Boltzmann distribution with the peak of the curve at the ground state.To fully characterize the ground state solution probability versus weight bit resolution, the results of running over 2000 random problems for weight bit resolutions from 1 to 5 is shown in Fig. 6(c).The ground state accuracy decreases monotonically as the bit resolution is increased.The total number of oscillator periods,   , required to reach ground state with 99% certainty is defined as   =   log(1 − 0.99) /log (1 − ()).
Here, () is the measured ground state probability from a single run for bit resolution b, and N c is the number of oscillator periods for a single run.We obtain the P(b) from Fig. 6(a) and plot the results in Fig. 6(d).Thus to achieve the ground state with 5 bits of weight resolution, an average of 13 periods would be required.
To characterize the cause of the non-ground state solutions for the non-binary weights, Fig. 7 plots the measured ground state probability versus the difference between the ground state Ising energy and the next highest energy state (Δ = 1 − 0) for randomly generated problems with 5 bit weights.Intuitively, this parameter characterizes the depth of the ground state potential well, with easier problems exhibiting large Δ, and more difficult problems exhibiting small Δ.The data shows that as Δ decreases, the probability of finding the ground state on a single trial also decreases.For example, when Δ < 2, the probability of achieving the ground state falls below 0.9.Noise in the system causes difficulty in distinguishing energynear solutions from the ground state solution.The minimum Δ value inversly scales with weight bit resolution, which is the cause of the decreased accuracy for higher bit problems.Alternatively, binary weight systems typically have very high ground state probabilities due to the large Δ values.For example, in the 4 node binary weight system, the minimum energy difference between E1 and E0, is Δ = 2, which gives rise to very high accuracy solutions.The complete absence of non-unitary probabilities for Δ > 2 implies that the system does not become trapped in local minimum, so long as the Δ is large enough.To compare the experimental results to simulations, we overlay a similar plot generated by simulating a set of randomly generated 5-bit weight problems with the SDE simulation technique described above.The simulated sample set includes 1000 different problems, and each point represents the mean ground state probability value of all the problems with a given Δ value.The behavior of the SDE simulated system matches reasonably well to the behavior measured experimentally.

Discussion
To understand the expected performance as the system size is scaled, the SDE solver discussed previously was extended beyond the case of V=4 and used to calculate the probability that the system can find the ground state as a function of V. We note that, for this scaling study, the annealing schedule described in Eq. ( 3) was used with τ = 5 cycles.This value was intentionally chosen to be small to allow for a pessimistic assumption of the solution time scaling trends.The problem set for this test was a set of circular graphs with each node connected to its two nearest neighbors and one node directly opposite (shown in Fig. 8(a)).
A deterministic solution to find the ground state of this problem set exists, which allows for the solution from the SDE solver to be compared against a global energy minimum.Fig. 8(a) plots the probability P that the system reaches the ground state solution after 500 independent runs for each problem set, showing a natural roll-off from ~100% ground state probability when V=8 to roughly 3% probability when V=300.
Translating these results to calculate the overall time required to find a solution for a given problem (  ), we first find Nt by defining Nc as the number of oscillation periods required for a single run of the simulated system to reach a minimum value and multiplying by the number of tries required to reach the ground state based on the expected probability P (as described in section 3c).Finally, multiplying Nt by the period of one oscillation yields the final solution time (  ),.Fig. 8(b) shows these computed values assuming a system built with 50 kHz oscillators (to match the current experimental circuit), along with a linear fit to a logarithmic trend line.The results are compared against the identical problem set and a similar analysis from previous Ising machine implementations 10 (orange points).For illustrative purposes, a similar trend line is plotted for a system assuming an oscillator frequency of 1 MHz.Direct scaling with oscillator frequency was demonstrated in Fig. 4, and the additional line is shown here to highlight the potential speedup that can be achieved against existing implementations with relatively modest hardware specifications.It is also important to note that this scaling analysis focuses entirely on ground state solutions.As shown in Fig. 6(b), when the ground state solution is not reached, the system often settles to a value very close to the absolute minimum, which can be sufficient for many optimization problems.

Conclusion
In summary we demonstrate the simulation, design, experiment, and characterization of a parallel, all-to-all connected coupled oscillator system which correctly maps to the weighted Ising model.The system is capable of successfully solving random MAX-CUT problems with 98% success probability with binary weights.The proof-of-concept breadboard experimental demonstration enabled a study of accuracy as a function of the interconnect weight bit-resolution.We show that the performance maps accurately to SDE simulations, and use these simulations to predict the system behavior at larger node sizes, which show favorable scaling when compared against existing optoelectronic for similar problem sets.
Recently, GPU implementations of a mean-field algorithm have shown a 20x speedup when compared against existing optoelectronic Ising machines 15 .The system presented here, if implemented with higher frequency oscillators, could prove competitive against these GPU implementations, and a full comparison against these instances for randomized and more densely connected problem sets will be the study of future work.
A significant advantage of the coupled oscillator system presented here lies in the fact that it can be scaled using existing low-cost hardware and lends itself especially well to an integrated circuit implementation.One major challenge, and an important area of future study, involves the interconnect architecture required to densely connect a large number of oscillators.While time-multiplexed approaches have been implemented in recent instances 9,10 , frequency-or code-division multiplexed systems could also provide interesting pathways to efficiently allocate connection resources.Another interesting area for study involves the use of measurement and feedback approaches to increase the ground state solution probability for larger problem sizes 24,25 .As transistor scaling begins to limit the progress in conventional digital computing architectures, the demand for alternative strategies continues to grow, and oscillator-based Ising machines could represent a novel and important platform to push the bounds of certain classes of computational problems.

Supplementary Material
All experiments have been performed using commercial off the shelf components, obtained from Digikey and assembled on solderless breadboards.Table 1 shows a detailed list of the circuit components, part numbers, and item description.

Fig. 1 .
Fig. 1.(a) Schematic of a fully connected 4-node system.(b) Solution states of the degenerate 4 node state.(c) Simulation of the oscillator phases.(d) Simulation statistical analysis of the six degenerate solution states.

Fig. 2 (
b).The coupling coefficient polarity is controlled by the polarity of the output of the differential summing amplifier to the nodes ViR and ViL.Digital potentiometers, R12-R34, are employed to control the individual gains of each of the input oscillator signals.The actual gain term is determined by the ratio of the feedback resistor, RFB, and the digital potentiometers.An Arduino microcontroller board was used to apply the digital I2C communication signals to the digital potentiometers.Precise tuning of the bias voltages for the oscillator and coupling circuit is critical to ensure accurate solution performance of the system.

Fig. 2 .
Fig. 2. (a) Circuit diagram of the LC oscillator circuit.The input coupling is inserted at the gates of transistors M3 and M4.(b) Multiply and accumulate cross-bar array composed of digital potentiometers.(c) Photograph of full breadboard system.The analog coupling coefficients from the Ising Hamiltonian (  ) are mapped linearly to the ratio of the gains between the various oscillators.The digital potentiometers employed here have 1024 tap points, with a maximum resistance of 20 kΩ.The conversion from the analytical coupling coefficients to the digital potentiometer values () is shown the following equation:

Fig. 3 .
Fig. 3. (a) System level diagram of the oscillator control system.Blue lines indicate digital signals while red lines indicate analog signals.(b) Schematic timing diagram of the oscillator phases during system operation.

Fig. 4 .
Fig. 4. Measured phase versus time data during the last weight change, J34, with (a) 50 kHz and (b) 5 kHz oscillators.A solution time of approximately 5 cycles is measured.The time-to-solution directly scales with oscillator frequency.

Fig. 5 .
Fig. 5. (a)-(e) Measured oscillator amplitude data at various binary weight configurations.The lines connecting the oscillators correspond to a weight of 1, while an empty space corresponds to a weight of 0. Time trace data of the oscillators under various binary weight configurations is shown in Fig. 5 to demonstrate the functioning system.The voltage amplitudes of each oscillator are shifted by 10 V to help distinguish the individual oscillators.All oscillation voltage amplitudes in this experiment range from 0 to

Fig. 6 .
Fig. 6.Experimentally measured probability of solution energy for (a) 1 bit and (b) 5 bit weight resolutions after more than 2000 trials with random weights.Solution energy of E0 represents the ground state solution.Error bars represent the standard deviation.(c) The measured ground state probability versus weight bit resolution.(d) Estimated solution time in number of oscillator periods based on a 99% confidence.The corresponding shaded regions represent the standard deviation.

Fig. 7 .
Fig. 7. Experimentally measured ground state probability versus solution energy difference between the ground state Ising energy and the next highest energy state (Δ = 1 − 0) for randomly generated 5 bit weights.

Fig. 8 .
Fig. 8. (a) Simulated ground state probability on a circular graph versus number of oscillators.(b) Scaling simulation of the coupled oscillator system with 50 kHz and 1MHz oscillator frequencies.For comparison, data extracted from previous Ising machine demonstrations for an identical problem 10 (orange dots) is plotted alongside the simulated data, illustrating a potential advantage of this coupled oscillator system.

Table 1 .
Circuit component information