A Note on Analyzing the Stability of Oscillator Ising Machines

The rich non-linear dynamics of the coupled oscillators (under second harmonic injection) can be leveraged to solve computationally hard problems in combinatorial optimization such as finding the ground state of the Ising Hamiltonian. While prior work on the stability of the so-called Oscillator Ising Machines (OIMs) has used the linearization method, in this letter, we present a complementary method to analyze stability using the second order derivative test of the energy / cost function. We establish the equivalence between the two methods, thus augmenting the tool kit for the design and implementation of OIMs.

Introduction: Ising Machines (OIMs) offer a markedly different paradigm for solving combinatorial optimizationa problem class that finds increasing importance in areas from machine learning (e.g., training of Boltzmann machines [1]) to applications in logistics such as route planning [2].Despite the ample use cases of such problems, solving them efficiently using digital computers and algorithms remains a challenge [3].Consequently, this has motivated the quest for alternate computing paradigms such as Ising Machines to accelerate the solution of such problems [4].
Networks of coupled oscillators offer a promising hardware platform to realize Ising Machines, commonly referred to as, Oscillator Ising Machines (OIMs) [5].Such a platform offers the advantage of being compact, scalable, low power, and compatible with commercial CMOS foundry process technology [6].The underlying design principle of such a platform exploits the natural equivalence between the global minima of the energy / cost function () describing a network of coupled oscillators under second harmonic injection and the global minima of the Where,  = ( 1 ,  2 ,  3 , … ,   ) represents the oscillator phases, and  and  s are the coupling and second harmonic injection strength, respectively.From equation ( 1) and ( 2), it can be deduced that, The dynamics described in equation ( 1) reveal that the system of equations has multiple fixed points, where ⅆ ⅆ = 0 (also, −(∇) = 0).In fact, every spin configuration ( 1 ,  2 ,  3 , … ,   ),   ∈ {0, } is a fixed point, resulting in 2  such fixed points.We note that there might be other fixed points defined by ( 1 ,  2 ,  3 , … ,   ) where   ∉ {0, }.The fixed points ( 1 ,  2 ,  3 , … ,   ) with the lowest energy represent the optimal solutions to the Ising Hamiltonian while the rest represent sub-optimal solutions.In our prior work [8], we analyzed the stability of the fixed points, corresponding to various spin configurations, of the OIM.Specifically, we showed that tuning the ratio of the coupling strength among the oscillators () and the strength of the second harmonic injection ( s ) has a dramatic impact on the stability of the globally optimal and sub-optimal fixed points.This analysis was performed by using the linearization method, i.e., by computing the Lyapunov exponents, which in this case, are defined by the eigenvalues of the Jacobian matrix ().A fixed point (globally optimal or not) is attractive if all eigenvalues at that point are negative; and unstable if at least one eigenvalue is positive.
Results: While the Jacobian analysis essentially entails working with the first order derivatives, the purpose of this work is to present an alternate approach, based on the second order derivates test of  (using the Hessian Matrix).We show that for a system whose dynamics are of the form −(∇)  = ⅆ  ⅆ for OIM), the stability of the fixed points can be analyzed using the eigenvalues of the Hessian Matrix (  ).A fixed point is attractive if all eigenvalues of   are positive.For the unstable case, if some or all of the eigenvalues are negative, then the fixed point is a saddle point, or a local maximum, respectively.To prove this, we establish the equivalence between the two methods.
The Hessian matrix for () can be computed as, Here, the symmetricity of the second derivatives ( Equation ( 5) reveals that in OIMs, the Jacobian matrix is half of the negative of the Hessian matrix ( = − 1 2   ).Consequently, the magnitude of the eigenvalues of the Jacobian matrix (  ) will be half of the magnitude of the eigenvalues of the Hessian matrix (   ) but with the opposite sign, i.e.,   = − 1 2    .This implies that the condition for a fixed point to be attractive, when analyzing the Hessian matrix, is that all eigenvalues be positive.When some eigenvalues are negative, it entails that the fixed point is a saddle point; when all eigenvalues are negative, it can be inferred that the fixed point is a maximum.For the general class of gradient descent systems defined by −(∇)  = ⅆ  ⅆ ( > 0),   = −   .

Conclusion:
In summary, we have presented a complementary method that utilizes the eigenvalues of the Hessian matrix to compute the stability of the fixed points in the OIM and is generally applicable to the broader class of gradient descent systems whose dynamics can be defined by -−(∇)  = ⅆ  ⅆ .Our approach shows that the stability of the fixed points in the OIM can be analyzed by identifying whether they represent a minimum of the energy function.
, equation (4) can now be expressed as: