Ising Machines based on networks of mutually coupled oscillators
Note, this is a work in progress, like a diary to formulate our thoughts and collect ideas, questions, hypothesis and results as we work towards gaining insight to the topic.
[collaborators: Dimitris Prousalis, Lucas Wetzel]
This post is about our exploration into the realm of so called Ising Machines. My curiosity into the topic was sparked by genius Bernd Ulmann, Mathematician, self-appointed Analog Evangelist, Inventor, Founder, museum director and most importantly a nice person.
Ising Machines are physical implementations of systems whose collective dynamics can be approximated by the Ising Hamiltonian. As shown, e.g., in “Ising formulations of many NP problems” by Andrew Lucas, many NP hard problems can be mapped onto the Ising Hamiltonian. Hence, if an Ising Machine can be build and scaled, it can act as an analog computer to solve these NP hard problems by finding/relaxing to the ground states of the system.
First of all I would like to link the first publications that I read to start with the topic. To be honest, they gave me a bit of a headache since by no means I was able to confirm the results that they presented when I stared looking at them in early 2021 — but more on that later. Here is the list to the open access manuscripts:
- OIM: Oscillator-based Ising Machines for Solving Combinatorial Optimisation Problems, Tianshi Wang, Jaijeet Roychowdhury
Hypotheses, issues and questions
- When scaling Ising Machines to large numbers of spins (represented by oscillators whose phases are binarized) and to operation at high frequencies, the systems dimensions may lead to non-negligible transmission time delays in the connections between the oscillators. We know, that such time delays can lead to multistability of in- and anti-phase synchronized states and therein. Furthermore, as the problem is usually encoded via the adjacency matrix of the network, these time delays will most likely be heterogeneous. Also, it has been shown that inert responses in, e.g., electronic oscillators, that stem from signal processing affect the self-organized dynamics in such systems.
The question is then, whether the mapping to the Ising Hamiltonian can still be justified and if not, which measures can be taken to achieve a sufficient (to solve a problem) approximation.
TODO: compile a list of open questions and problems related to the presence of time delayed coupling and inertia, quantify their effects in simulations and look for ways to address these.
- When simulating delay-coupled phase oscillators using the Kuramoto model we made the following observations:
a) The results in Oscillator-based Ising Machine, chapter IV: Examples can be reproduced for sufficiently strong dynamical noise. If the noise was too small, we did not observe the convergence to the correct MAX-CUT solutions in all cases (actually we observed that only rarely). This could be a sign for rugged energy landscapes that require some ‘shaking’ to find the absolute minima and hence the correct solutions to the MAX-CUT problem.
b) When simulating using a second order Kuramoto model that takes into account inertia like effects introduced by, e.g., a loop filter in a phase-locked loop or the mass of a pendulum, the phases did not converge to the correct solutions for arbitrary loop filter cut-off frequencies. Increasing the cut-off frequency, or in other words, reducing the oscillators inert behavior, led again to the correct solutions. My hypothesis is, that loop filters with cut-off frequencies far below the second harmonic injection frequency will damped this too strongly and prevent he binarization of the oscillators phases. Hence, it may be beneficial to mix the second harmonic injection locking signal directly into the control signal of the PLL. In principle it may be possible to use a ‘copy’ of the feedback signal which is mixed with itself. The resulting output can then be filtered by a band-pass that allows only the second harmonic. These are however details the implementation of such an idea, in e.g., a PLL type of setup.
c) Adding small transmission delays seemed not to hinder convergence to the correct solutions, independently of whether the time delay would have led to in- or anti-phase synchronized states.
Playground, simulations & observations
Here I collect some results obtained from simulating different types of Kuramoto models. The code used here can be found on Github. The documentation process is still underway and I am being supported in cleaning up and organizing the code by Deborah Schmidt.
Along the results presented in Oscillator-based Ising Machine obtained by Matlab simulations of the Kuramoto model we start here and ask whether these results can be reproduced. Our model reads in the most general form [1,2,3]:
where denotes the instantaneous frequency, and the intrinsic frequency of oscillator , with . denotes the coupling strength or sensitivity towards input signals to oscillator from other oscillators in the network, as given by , the components of the adjacency matrix that are either one or zero depending on whether there is a connection from oscillator to or not, respectively. The integral represents the signal processing or filtering within the feed-forward path of the oscillators, e.g., a low pass filter in a phase-locked loop oscillator with impulse response function . denotes the phase of oscillator , the transmission time delay for signals propagating from oscillator to , and the feedback time delay within oscillator . denotes the strength of the harmonic injection locking. Finally, denotes additive Gaussian white noise with variance .
In terms of the model, we aim to address how the following parameters affect finding the ground states and relate to the classical Ising model:
- how to appropriately increase or , which one preferably and at what rate, does it depend on the network size?
- the role of the loop filter via the , the filter transfer function – if the cut-off frequency is below the harmonic frequency, then the injection locking signal will be damped/filtered and binarization may be impaired (here we should consider a modified PLL architecture where the harmonic injected signal does not pass trough the loop filter, e.g., multiplication by itself and band-pass to extract harmonic)
- the role of feedback and transmission time delays, as the size of a system grows it becomes increasingly difficult to arrange the oscillators such that they can be coupled to any other in the network with the same transmission time delay; question: if there are transmission time delays that are heterogeneous within one connection (asymmetric), i.e., between an oscillator and and oscillator , this leads to phase-difference between these oscillators for the case of pure mutual coupling without injection locking – how does that affect a networks ability to relax to a ground state of the Ising model? the same holds for the case and
- what about the coupling capacity/normalization of the inputs? does an oscillator couple more strongly the more inputs it receives?
- should the coupling between the oscillators be repelling, i.e., should ?
In the publication Oscillator-based Ising Machine, see chapter IV, they set , and simulate different network topologies, using simple MAX-CUT problems to benchmark. Note, that the coupling strength is increased linearly / ramped up once the simulation has started.
Using the same parameters we find the following results for a system of oscillators, where oscillator receives input from all other oscillators while all other only receive from oscillator . Hence, we have a star-like topology with oscillator being the center. In this case, the MAX-CUT problem is solved if all oscillators evolve to a phase different by from that of oscillator . The results were promising and lead to the correct solution of the MAX-CUT problem in most cases, however not in all cases as we show in the figures below:
We then realized the importance of the process of ramping up the coupling strength for finding the correct solution successfully. Also, adding dynamical frequency noise seemed to help convergence. Moreover, we simulated some realizations with time delay, signal filtering and ramping up the harmonic injection coupling strength instead of the coupling strength of the first harmonic coupling terms. This had profound impact on whether the correct solution to the optimization problem was found. Therefore we decided to develop ideas on how to approach the impact of these parameters in a more structured way. Our preliminary results is what we will present in the following.
How we measure the success probability and time to solution as we change different system parameters?
We decided to run a number of realization for a fixed set of parameters and solve the MAX-CUT problem for a system of oscillators arranged in a star like topology, where an oscillator indexed by is located in the middle of the star and hence receives input from all others, while all others only receive input from , the center oscillator. After computing realizations we ask how many have led to a correct solution and in what average time that was achieved. A realization is categorized at being solved correctly when the magnitude of the order parameter settles at a value of in that particular case, i.e., all phases have separated correctly. The Kuramoto order parameter is defined as
where denotes the magnitude of the order parameter and the mean phase. measures the synchrony of the oscillators and approaches if all oscillators are in-phase. In our benchmark case, all outer oscillators of the star topology form a subset that is shifted by a phase difference of with respect to oscillator . One can then calculate
which yields the result for the particular MAX-CUT problem used for benchmarking here mentioned above and visible in Fig.~1.
A few more words on how we obtained the statistics on the realizations automatically from the simulations.
In order to determine whether the correct solution to MAX-CUT problem was found we checked the asymptotic values of the phase differences. Hence, if the phase differences of all oscillators with respect to the phase of are equal to asymptotically, the correct solution has been obtained. We quantify this by checking these phase-differences and the value of the order parameter. Since the phase-configuration should be constant once the system has settled to the ground state, we compute the average over one and a half period of the intrinsic oscillations of the oscillators.
Extracting the mean time until a solution is obtained is measured as follows. We consider only those realizations that led to a correct solution of the MAX-CUT problem. For these, we calculate the time-derivative of the order parameter which we then average over half a period of the intrinsic frequency of the oscillators. This yields a more smooth curve. We then track from the end of the time-series towards the beginning, i.e., backwards in time, when the order parameter changes for the first time. This represents the end of the transient dynamics that led to the solution. To calculate the time to solution, we choose the difference between the time when the transient dynamics have finished and the time where we started to ramp up the coupling strength.
Here is an example how the resulting plots look like: