Designing Ising machines with higher order spin interactions and their application in solving combinatorial optimization

The Ising model provides a natural mapping for many computationally hard combinatorial optimization problems (COPs). Consequently, dynamical system-inspired computing models and hardware platforms that minimize the Ising Hamiltonian, have recently been proposed as a potential candidate for solving COPs, with the promise of significant performance benefit. However, prior work on designing dynamical systems as Ising machines has primarily considered quadratic interactions among the nodes. Dynamical systems and models considering higher order interactions among the Ising spins remain largely unexplored, particularly for applications in computing. Therefore, in this work, we propose Ising spin-based dynamical systems that consider higher order (> 2) interactions among the Ising spins, which subsequently, enables us to develop computational models to directly solve many COPs that entail such higher order interactions (i.e., COPs on hypergraphs). Specifically, we demonstrate our approach by developing dynamical systems to compute the solution for the Boolean NAE-K-SAT (K ≥ 4) problem as well as solve the Max-K-Cut of a hypergraph. Our work advances the potential of the physics-inspired ‘toolbox’ for solving COPs.

The minimization of the Ising Hamiltonian using dynamical systems such as coupled electronic [1][2][3][4][5] and photonic oscillators [6][7][8] has received substantial attention in recent years 9,10 . A significant driving force behind the effort to realize a so-called 'Ising machine' is that the solution to the Ising model can be mapped to many computationally intractable problems in combinatorial optimization (e.g., MaxCut, Traveling Salesman Problem (TSP) among others) [11][12][13][14][15][16][17][18] . Consequently, this creates the possibility of realizing Ising machine-inspired custom accelerators that can offer the possibility of significant performance benefits. However, dynamical system formulations that have been used to 'solve' the Ising model typically consider only pair-wise coupling; examples include, oscillator Ising machines, coherent Ising machines etc. From an application standpoint, while these characteristics capture quadratic interactions, the dynamical systems and their supporting computational models cannot be applied directly to solve problems that require higher order interaction among the spins 19,20 . Therefore, the objective of this work is two-fold: (1) define dynamical systems that model higher order (> 2) interactions among the Ising spins; and (2) map the resulting dynamics to relevant computational problems. We consider two examples: computing the solutions for the NAE-K-SAT (Not-All-Equal SAT) problem and the Max-K-Cut of a hypergraph. Our motivation behind selecting these two combinatorial optimization problems was that their objective functions directly map to the solution of the higher order Ising models, and therefore, help illustrate the principle of how dynamical systems for the higher order Ising models can be used in combinatorial optimization. Also, we emphasize here that presently our focus is on defining the Ising machine dynamics that capture the higher order interactions, and not on the physical implementation of the higher order interactions.
The general form to represent higher order interactions among the Ising spins can be expressed as, where J ij represents the pairwise interaction coefficient between two Ising spins. The first term on the right-hand side ( − i,j J (2) ij s i s j ) is usually considered when describing quadratic/pairwise interactions among Ising spins www.nature.com/scientificreports/ s = {−1, 1} n ; the Zeeman term which considers the interaction of spins with an external magnetic field has been neglected here. Considering the higher order interactions among the spins can help describe the objective functions of several combinatorial optimization problems (COPs) as illustrated here with the example of the NAE-K-SAT problem (without the need for problem decomposition). The NAE-K-SAT problem is a constrained version of the Boolean Satisfiability (SAT) problem where the objective is to find an assignment for the variables of the given Boolean expression (in the conjunctive normal form) such that: (a) at least one variable in every clause is TRUE (i.e., the clause is satisfied; standard SAT constraint); (b) at least one variable in every clause is FALSE 21 ; the NAE-K-SAT problem is considered here since it directly maps to the general form of Eq. (1), as illustrated further on. Using an approach inspired by SAT, the NAE-K-SAT problem can be expressed as computing an assignment for the variables such that Y (= C 1 . . . x N ) (i.e., S i and C i have the same variables but in opposite forms). Traditionally, when considering only pairwise interactions among the Ising spins, mapping such problems can entail significant pre-processing including the use of auxiliary variables that can significantly increase the size of the problem that must be eventually solved 20,22-25 using the dynamical system.

NAE-4-SAT.
To illustrate how we can map the NAE-K-SAT problem to higher order interactions among the Ising spins, we first consider the example of the NAE-4-SAT problem where each clause of the NAE-4-SAT problem consists of 4 literals, expressed in the general form as x is a set of Boolean variables). K = 4 is specifically chosen here since it is the lowest K where higher order interactions among the Ising spins are required to formulate the objective function for the problem (shown in Table 1). To formulate the problem in terms of Ising spins, we utilize the following property among the Boolean variables and the spins Here, the logic level 0 (1) corresponds to an evaluation of − 1(1) of the expression on the right-hand side, respectively. Furthermore, the complement of the logical OR among the XOR terms ( can be expressed as, 1 − . . . 1+s k s l 2 ≡ 1 8 (1 + s i s j + s i s k + s i s l + s j s k +s j s l + s k s l + s i s j s k s l ) . It can be observed that besides the second order interaction terms, the resulting expression also contains a 4 th order interaction term among the spins. Consequently, the objective function for the NAE-4-SAT problem, over M clauses, can be formulated as the minimization of Table 1. Objective functions for the NAE-K-SAT problem expressed using Ising spins. We note that constants and scalars have not been shown here in the expression for the single clause as well as for the objective function.

2
Expression for a single clause: It can be observed that when the variables appear only in the normal form i.e., c mi ≥ 0 , the expression represents the solution to the archetypal MaxCut problem 3 Expression for a single clause: Expression for a single clause: Objective function: Expression for a single clause: ≡ s i s j + s i s k + s i s l + s i s m + s j s k + s j s l + s j s m + s k s l + s k s m + s l s m + s i s j s k s l + s i s j s k s m + s i s j s l s m + s i s k s l s m + s j s k s l s m Objective function: www.nature.com/scientificreports/ Here, c mi = 1(−1) , if the ith variable appears in the mth clause in the normal (negated) form; c mi = 0 if the ith variable is absent from the mth clause. Using the same approach, we derive such expressions for a few other values of K in the NAE-K-SAT problem in Table 1. Details of the derivation of the objective function for NAE-5-SAT are shown in Supplementary 1. Constructing a dynamical system for the NAE-K-SAT problem. We now aim to formulate the dynamical system and the corresponding energy function for the NAE-K-SAT problem. The dynamical system, defined by − ∇ φ E i = dφ i dt , is designed such that the ground state of the 'energy' function (more precisely, the Lyapunov function) must correspond to a global optimum of the objective function. To construct this system, we draw inspiration from the dynamics of coupled oscillators under second harmonic injection which effectively forces the oscillator states to assume a binary phase value of 0 or π (details of the second harmonic injection can be found in work by Wang et al. 26 ). Without loss of generality, we assume that one spin state (say, s = +1 ) is represented by phase 0 while the other spin state ( s = −1 ) is represented by the phase angle π. Subsequently, the second order interaction terms among the Ising spins s i s j can be represented by cos(φ i − φ j ) . When the spins are in opposite states i.e., s i = 1(−1); s j = −1(1) , s j s j ≡ cos φ i − φ j = −1 , whereas when the spins are in the same states i.e., s i = 1(−1); s j = 1(−1) , s j s j ≡ cos φ i − φ j = 1 . Similarly, the higher order interactions can be modeled as shown in Table 2.
The equivalence between the higher order terms and the corresponding energy term is shown in Table 3.
Using the above relationships developed in Table 1, the energy functions for the NAE-K-SAT problem can be formulated as shown in Table 4. The corresponding dynamics dφ i dt , shown in Table 4, can be obtained from the dynamical system equation is added to ensure that the oscillator phases effectively binarize to {0, π} . The energy contribution of this term is minimized ( = −N C s 2 ) at the binary phase points φ ∈ {0, π} . Consequently, by using the appropriate strength of the second harmonic injection ( C s ), we can ensure that the energy function reaches its minimum for φ ∈ {0, π} . We have borrowed this approach from prior work on oscillator-based Ising machines (with second order interactions) 26 .
Furthermore, using the dynamical system equation dφ i dt = − ∇ φ E i , we can also show that for the energy functions described in Table 4, dE dt ≤ 0 i.e., they are Lyapunov functions. Figure 1 shows an illustrative example of the NAE-4-SAT problem computed using the proposed dynamical system. Details of the simulation used to simulate the illustrative NAE-4-SAT problem are described in Supplementary 4.
Max-K-Cut on a hypergraph. In the prior section, we exploited the binary nature of the Ising spins (along with higher order interactions among them). We now 'extend' the definition of the 'spin' in order to facilitate the design of computational models for an even broader spectrum of COPs that would benefit from the use of > 2 states for each node/spin. To facilitate this, we express the possible states of a spin as re iθ k , where r = 1 , and θ k = 2πk K ; k = 1, 2, . . . K − 1 . When K = 2 , the possible states are within {1, − 1}, which represents the traditional definition of an Ising spin. In contrast, when K > 2 , the 'spin' assumes K configurations, represented as complex quantities (e.g., for K = 3 , the possible states are 1 , e i 2π (1) 3 , e i 2π (2) 3 ). While we have utilized this concept for solving combinatorial problems on graphs (i.e., problems with quadratic objective functions) 18 , here we  www.nature.com/scientificreports/ explore this concept for hypergraphs (that entail higher order interactions) by considering the example of solving the Max-K-Cut of a hypergraph. Computing the Max-K-Cut on a hypergraph is defined as the challenge of partitioning the nodes of a hypergraph into K partitions in a manner that maximizes the number of hyperedges having nodes that lie in at least two sets created by the partitions 27 . The Max-K-Cut problem and its comparison with the archetypal MaxCut problem are illustrated in Fig. 2a,b for the case of a graph and a hypergraph, respectively.
To develop the objective function for the problem, each hyperedge of the graph can be expressed as , where s j = 1e iθ j ; θ j can assume any of the following values from 2πk K ; k = 1, 2, . . . K − 1 enforced by the higher order harmonic injection. c mj = 1(0) if the jth node belongs (does not belong) to the mth hyperedge. We note that the ′ i ′ represents the imaginary number √ −1 whereas ′ i ′ refers to the index.
f (�θ ij ) is designed such that Re s i s * j e if (�θ ij ) = −1(1) , if the nodes i and j are placed in different (same) sets, and essentially rewards (penalizes) the system in terms of energy, respectively. Additional details about the design and properties of f (�θ ij ) have been presented in our prior work 18 . Consequently, if the hyperedge satisfies the Table 3. Equivalence between the higher order Ising spin interaction terms and the equivalent energy function.
Second order interactions ( s i .s j ) Third order interactions ( s i .s j .s k ) Fourth order interactions ( s i .s j .s k .s l ) www.nature.com/scientificreports/ Table 4. Objective functions, corresponding energy expressions, and system dynamics for NAE-K-SAT problems for K = 2, 3, 4, and 5. We note that while the form of the expressions for K = 2 and K = 3, as well as K = 4 and K = 5 are similar, the coefficients ( c mi ) are different. C is the strength of coupling among the nodes whereas C s represents the strength of the second harmonic injection. Objective function: Energy function: www.nature.com/scientificreports/ criterion for the Max-K-Cut i.e., that the nodes that are connected by it belong to at least two sets, the corresponding h m assumes a value of 0 , else h m = 1 . Subsequently, the objective function for the problem, which entails maximizing the number of such hyperedges, can be expressed as minimizing H , where, www.nature.com/scientificreports/ As an example, considering a hypergraph where the maximum number of nodes connected by a hyperedge is 3, the objective function for the Max-K-Cut problem can be expressed as: where, For a hypergraph with hyperedges having more than 3 nodes, the objective function entails the use of higher order interactions among the spins.
To formulate a dynamical system for minimizing the above objective function, we express Re s i s * j e if (�θij) as cos �θ ij + f �θ ij . Furthermore, we restrict the configuration space of θ to 2πk K where k = 1, 2, . . . K − 1 , by injecting the Kth harmonic (of sufficient strength) which lowers the energy at specific phase points, as described in prior work 18 . The resulting energy function can be described as, We note that φ has been used to express the energy function for the dynamical system instead of θ which represents the configuration space of the 'extended spin' . The corresponding dynamics for which the function in Eq. (8) is a Lyapunov function are given by: In the derivation of Eq. (9b), we exploit the fact that ∂f (�φij) ∂φ i = 0 18 . Furthermore, using Eq. (9a), it can be ≤ 0 (similar to Eq. (3)). We now evaluate our proposed model on a representative hypergraph. We consider a hypergraph where each hyperedge has 3 vertices. The corresponding dynamics for this case can then be written as, Figure 2c-k also shows the computed Max-K-Cut (for K = 2, 3, and 4) for a hypergraph instance (with 10 nodes, and 20 hyperedges). The illustrative problem has a maximum of 4 nodes per hyperedge. Details of the simulation used to simulate the illustrative Max-K-Cut problem are described in Supplementary 4.

Conclusion
In this work, we develop computational models for Ising machines that consider higher order interactions (beyond quadratic/pairwise) among Ising spins. Our approach enables the direct formulation of analog computing models for many COPs that entail such interactions without the need for problem decomposition and reduction. Furthermore, using the combination of higher order interactions along with 'expanding' the number of 'spin' states to greater than 2, we can directly map and solve an even broader class of problems on hypergraphs. This has been summarized in Fig. 3. While the focus of the present work was to develop dynamical systems as higher order Ising machines, evaluating the scalability of this approach i.e., its ability to solve larger graphs will be crucial to its eventual success and utility. As the graph sizes increase, the role of local minima in the high dimensional phase space becomes increasingly important. The system dynamics may get trapped in such minima resulting in sub-optimal solutions. Furthermore, identifying the optimal range for the parameters ( C and C s ) in larger systems may also become more challenging. Eventually, these factors will also ascertain the performance benefits of this approach over traditional digital algorithms used to solve such problems. A systematic study to evaluate the scalability of this approach and its comparison with digital methods will be undertaken in the future.
In the context of the broader effort focused on developing dynamical system-inspired models for solving hard COPs, this work expands on the potential of physics-inspired solvers to accelerate COPs.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Received: 24 January 2023; Accepted: 5 June 2023 Figure 3. Proposed classification of COPs based on the number of states/configurations for the nodes, and the nature of interaction among them. Categorizing COPs helps develop a framework to formulate dynamical systems to solve the COPs. The present work enables the direct development of computational models for COPs that entail higher order interactions.