Efficient optimization with higher-order Ising machines

A prominent approach to solving combinatorial optimization problems on parallel hardware is Ising machines, i.e., hardware implementations of networks of interacting binary spin variables. Most Ising machines leverage second-order interactions although important classes of optimization problems, such as satisfiability problems, map more seamlessly to Ising networks with higher-order interactions. Here, we demonstrate that higher-order Ising machines can solve satisfiability problems more resource-efficiently in terms of the number of spin variables and their connections when compared to traditional second-order Ising machines. Further, our results show on a benchmark dataset of Boolean k-satisfiability problems that higher-order Ising machines implemented with coupled oscillators rapidly find solutions that are better than second-order Ising machines, thus, improving the current state-of-the-art for Ising machines.


Introduction
An Ising machine is a type of parallel computer utilizing energy relaxation in a network of interacting binary variables. Ising machines have been proposed as efficient methods for finding optimal or near-optimal solutions to hard combinatorial optimization problems [1]- [6]. For a given combinatorial optimization problem, the network interactions are shaped so that the energy minima correspond to the problem solutions. For mapping a given combinatorial optimization problem to a network, a common strategy is to formulate the objective as the energy function of an Ising model, an abstract network of coupled bipolar variables originally proposed to model ferromagnetic material. The Ising model can then be implemented on hardware, referred to as an Ising machine. Ising machines implemented on quantum computers promise optimal solutions [3], [7]- [9]. However, due to the challenges of constructing them, Ising machines based on classical physics are reemerging and new technologies are being developed. There is a large variety of possibilities for implementing classical Ising machines, including coupled electrical oscillators [6], [10]- [12], optical parametric oscillators [13], stochastic circuits (probabilistic bits) [14], and neuromorphic hardware [1], [15], [16]. Here, we focus on classical Ising machines for approximately solving combinatorial optimization problems at scale and extremely fast.
Casting a combinatorial optimization problem as an Ising model usually takes two or three steps. The first step is to express the combinatorial optimization problem objective as a polynomial in the binary variables. The second step is mapping the polynomial to the energy function of an Ising model. For many combinatorial optimization problems, step one results in a higher-order polynomial [2], [17]- [21], i.e., a polynomial with terms that contain products of more than two binary variables. However, most Ising machines utilize second-order polynomial interactions between variables. In this case, a third step, called quadratization [2], [17], [18], [20], [22]- [25], is applied for reducing higher-order terms in the polynomial to second-order. The resulting second-order polynomial represents the energy function of a classical Ising model, i.e, a second-order network in which each interaction just couples a pair of variables [26]. Quadratization increases the network size by adding auxiliary variables and it requires increased precision and range of the second-order interaction coefficients compared to higher-order interactions [18], [19].
Higher-order Ising models -models that include polynomial interactions of a degree greater than two -have received little attention because the possible number of interactions grows exponentially with the interaction order. Thus, the training and implementation of higher-order Ising models seemed intractable and impractical [27]. Here, we propose to skip the step of quadratization and instead use higherorder Ising models that directly implement the higher-order polynomials describing the combinatorial optimization problems. Although this proposal seems daunting at first glance, we show that for important classes of combinatorial optimization problems, the corresponding higher-order Ising machines require fewer variables and connections than the second-order Ising machines resulting from the quadratization approach.
Among the proposed Ising machines, coupled electrical oscillators are promising for combinatorial optimization problems in terms of solution quality [28], and the ability to leverage existing technologies such as complementary metal-oxide-semiconductor (CMOS) ring oscillators [29], [30]. To build an oscillator Ising machine, the continuous phases of oscillator variables have to be biased towards two anti-symmetric states, for example, by sub-harmonic injection locking [28], [31], [32]. To demonstrate a concrete higherorder Ising machine, we investigate a network of coupled Hopf oscillators with sub-harmonic injection locking, referred to as a higher-order oscillator Ising machine. Results from our simulations show that the higher-order oscillator Ising machine not only uses fewer network resources compared to the second-order oscillator Ising machine but, importantly, achieves better solutions. All told, our results suggest that, against common beliefs, optimization with higher-order Ising machines can outperform traditional Ising model approaches.

Mapping constraint satisfaction problems to Ising models
A broad class of combinatorial optimization problems are constraint satisfaction problems, including invertible logic circuits, Boolean satisfiability (SAT) problems, and Boolean maximum satisfiability (MaxSAT) problems. SAT solvers have many direct applications in areas, such as artificial intelligence [33], electronic design automation [34], cryptography [35], and many more. Many Boolean constraint satisfaction problems naturally map to higher-order polynomials [2], [17]. The most common approach for solving constraint satisfaction problems with Ising machines has been first to apply quadratization for translating problems to second-order polynomials, and then use second-order Ising machines to solve them efficiently [2], [18], [19], [22], [24], [26]. However, optimization can also be performed in higher-order Ising machines without quadratization [21], [36], [37]. Here, we aim to construct higherorder Ising machines for Boolean constraint satisfaction problems which are simple, yet, scale to large problems and quickly find near-optimal solutions.
In Boolean constraint satisfaction problems, the Boolean variables must take a state which satisfies a set of pre-defined constraints. For the hth constraint containing k variables, the state space, S h = {−1, 1} k , can be partitioned into two sets. Let C h = {c ∈ S h : c = satisfied state} be the set of valid states, i.e., that satisfy the constraint, andC h = S h \ C h = {c ∈ S h : c = unsatisfied state} be the set of invalid states which do not satisfy the constraint. Any logic function can be expressed by a constraint for which the set C h represents the truth table of the function. An objective or energy function of the hth constraint, E h , can be written as the characteristic function of its set of invalid states [2]: or, equivalently (Methods 4.2), as one minus the characteristic function of its set of valid states: Thus, the sizes of the sets of valid and invalid states may determine which of the two equations is preferable. Let N C h = |C h | and NC h = |C h | denote the size of the set C h andC h , respectively. Then, Eqs. (1) and (2) contain a sum with NC h and N C h terms, respectively. Note that both energies contain higher-order interactions of the order of the size of the constraint.
The total energy for a constraint satisfaction problem is the weighted sum of the individual constraints, Eq. (3): Eq. (3) generalizes our method to weighted MaxSAT problems, which have many applications [38]. In MaxSAT, each constraint is assigned a weight, w h , representing the relative importance of satisfying the hth constraint. Here, Γ is the set of indices for the problem constraints, E h is energy function for the hth constraint formulated according to either Eq. (1) or (2).
b d e f g h c Figure 1: Mapping optimization problems to Ising models. Two example problems. Left: XOR circuit. a, Circuit schematic for the XOR. The XOR gate has two inputs, s 1 and s 2 , and one output, s 3 . b, State table has eight lines. Four lines are input configurations for valid/true output, the other four lines input configurations for invalid/false output. c, The higher-order energy function for the XOR in both the factored and simplified form. d, Energy and corresponding hypergraph of third-order XOR Ising network, variables nodes, depicted as circles, connected by one interaction, depicted as a square. e, Energy and corresponding graph of second-order XOR Ising network, resulting from quadratization. The graph contains four variable nodes (one auxiliary variable), six second-order interactions and four first-order interactions (biases). Right: SAT problem. f, SAT problem in CNF. The SAT function is written with binary variables, x i ∈ {0, 1}, wherex i denotes the variable negation. g, The SAT problem in CNF has an equivalent circuit representation consisting of k-input OR gates which output to one AND gate. h, The energy can be succinctly formulated with one term per clause using Eq. (1).
Eq. (1) or (2) are higher-order interactions represented as factored polynomials. The equations can be expanded to coincide with the common formulation of a higher-order Ising model (Eq. 5 in Methods 4.1). Either the factored or expanded parameterization may be preferred depending on the problem and which form results in the fewest number of terms in the energy. In general, the expanded energy may contain 2 k − 1 terms or parameters. However, for many practical problems each clause contains only a few literals, hence, k is small. The factored representations require N C h k and NC h k parameters for Eqs. (2) and (1), respectively. Thus, when k is large or the expanded form does not simplify to a few terms, the factored representation is preferable.
The derivation of Ising models is first explained for two small examples of combinatorial optimization problems, the exclusive OR (XOR) invertible logic gate, and a small SAT problem. The XOR problem can be depicted by the XOR gate symbol (Fig. 1a), and its state table (Fig. 1b). The expanded and simplified energy polynomial of XOR contains only one interaction (Fig. 1c), resulting in a very simple hypergraph of the corresponding third-order Ising network (Fig. 1d). The quadratization of the thirdorder XOR polynomial produces a second-order Ising network with one additional auxiliary variable, six second-order interactions, and four biases (Fig. 1e). The additional network resources required after quadratization may be negligible for small problems but significantly change the scaling behavior of required resources for larger problems (Fig. 2).
Any SAT problem can be written as the product (conjunction or AND) of clauses (constraints) where each clause is the Boolean sum (OR) of literals. A literal is a variable or its negation. This form is known as conjunctive normal form (CNF). For a particular 3 clause SAT problem, the CNF ( Fig. 1f ) corresponds to a logic gate circuit (Fig. 1g), and a factored higher-order energy polynomial (Fig. 1h). The factored energy polynomial of a SAT problem corresponds to Eq. 3 with w h = 1 ∀h. Therefore, any SAT problem in CNF maps directly to a higher-order Ising model in which each higher-order interaction represents a clause. The order of an interaction corresponds to the size of the corresponding clause.

Model scaling of higher-order and traditional Ising models
Quadratization of higher-order interactions introduces auxiliary variables and adds second-order interactions (XOR example in Fig. 1), thereby potentially increasing the total resources required by the corresponding Ising machine. To quantify this effect, Fig. 2 compares the resource use of higher-order models versus second-order models on kSAT benchmarks [39], [40]. k SAT is a SAT problem where each clause involves maximally k variables. Quadratization of kSAT proceeds first by reducing a kSAT problem to 3SAT for k > 3, which can always be done [41], and then quadratization of the 3SAT problem. We use the D-Wave Ocean software package for quadratization (Methods 4.5), which accepts the minimum classical energy gap, ∆E min , as an input parameter. ∆E min is the difference in energy between satisfied states and the lowest energy unsatisfied state. The choice of minimum energy gap value influences the annealing time in quantum adiabatic annealing [3] and the state acceptance probability in simulated annealing [42]. Increasing the minimum energy gap for an Ising machine may improve the optimization, however, it tends to increase the number of auxiliary variables and interactions required (Methods 4.5). We compare higher-order to second-order models in terms of the number of variables in the energy function and the number of connections needed to implement all interactions. We consider second-order models with different minimum energy gap values. Nearby values of ∆E min result in the same quadratization, therefore, we investigate ∆E min settings of 1, 5, 10, and 13 where 1, 2, 3, and 5 auxiliary variables are introduced per clause, respectively. In addition, we found that the method used to perform quadratization increases the required precision or resolution of coupling coefficients from one bit for factored higher-order Ising models to at most six bits. This is another significant difference in resource requirements, as hardware typically offer limited resolution precision for representing interactions [30].
To compare the resource use of interactions of different orders we consider the number of connections between nodes that are required for their implementation. The required number of connections depends on the way a higher-order interaction is implemented, here we compare two methods of implementation. The first method is bidirectional connections between all variables participating in the higher-order interaction -a k th-order interaction requires k(k − 1) connections (Fig. 2c). The second method uses an intermediate computational node that receives input from all other variables participating in the interaction and sends output back to all other variables -a k th-order interaction requires 2k connections (Fig. 2e).
Our comparison shows that second-order models based on quadratization of higher-order models require a much greater number of variables and connections compared to higher-order models for kSAT benchmarks Fig. 2. In particular, second-order models require three orders of magnitude more variables and one order of magnitude more connections compared to higher-order models. In addition, the number of variables obtained from the D-Wave Ocean software package for ∆E min = 1 is the same as another method of quadratization based on a circuit decomposition of SAT clauses [43] which introduces one auxiliary variable per clause (Methods 4.5).

Solving SAT problems with a higher-order oscillator Ising machine
To compare the computation performance of higher-order Ising machines with corresponding secondorder models, we use a concrete network model of coupled oscillators. In our higher-order oscillator Ising machine, each oscillator is described by a complex variable z i ∈ C, which evolves according to: Here, z i represents the amplitude and phase of the ith oscillator. On the right-hand side, f (z i ) is the local oscillator dynamics, and ∂E(g(z(t))) ∂zi the partial derivative of the Ising energy with respect to oscillator z i , with time-dependent coupling coefficient r i (t), and optional element-wise non-linearity, g(z(t)) = z(t)/|z(t)| for normalizing the amplitude of each oscillator. Further, l(z i ) =z i is the phase quantization  The bars are grouped by k, the number of variables per clause. a, The ratio of the number of variables required for second-order networks compared to higher-order networks. c, A higher-order interaction implemented with all-to-all connectivity. e, A higher-order interaction is implemented with a computational node for each constraint. b, d, The ratio of the number of connections required for second-order networks compared to higher-order networks implemented with all-to-all connectivity, c, and intermediate computational nodes, e.
signal driving the phase of oscillator z i to discrete states, with time-dependent "annealing" coefficient, q i (t). The phase quantization signal is equivalent to sub-harmonic injection locking (Methods 4.7).
We compared simulations of higher-order oscillator Ising machines and second-order oscillator Ising machines for solving kSAT benchmarks [39], [40]. Higher-order oscillator Ising machines achieve better solutions than second-order oscillator Ising machines on all 3SAT benchmark problems, as measured by mean energy at the solution points (Fig. 3a). Only for the smallest problem instances (20 variables), the difference is small. For larger problems, a substantial gap in energy appears and increases with problem size. Interestingly, even second-order oscillator Ising machines with large minimum energy gaps and, correspondingly, high resource use cannot close the performance gap to higher-order oscillator Ising machines. The performance gap amounts to about 0.75 percent of constraints satisfied for the large 3SAT problems (Fig. 3b). Finding optimal solutions, i.e, states which satisfy all the constraints, is a hard problem as there could be very few satisfying states in the entire state space. Nevertheless, for larger problems of the 3SAT benchmarks, higher-order oscillator Ising machines tend to find solutions that satisfy all constraints with greater probability than the second-order oscillator Ising machines, Fig. 3c. In fact, the higher-order oscillator Ising machine is the first reported Ising machine to find satisfiable solutions to the largest 3SAT problems (250 variables) since the previous efforts with second-order Ising machines have been unable to find solutions satisfying all clauses [43], note the missing bars in Fig. 3c.
Annealing typically improves the quality of solutions found by Ising machines [20], [21], [28], [43]. In both our higher-order oscillator Ising machine and existing second-order oscillator Ising machines, a process analogous to adiabatic and simulated annealing is achieved by gradually increasing the coefficient in the sub-harmonic injection locking term, q i [10]. We investigated linear annealing schedules with different duration, measured by the number of cycles of the resonant frequencies of the oscillators. The percentage of constraints satisfied at the end of the annealing schedule improves with the duration of the annealing schedule (Fig. 3d). The time-to-solution for reaching a fixed target of 95% of constraints satisfied (TTS 95 ) scales linearly with the slope in the annealing schedule (Fig. 3e). For large slopes, the TTS 95 can be a fraction of a cycle, consistent with previous findings that oscillator Ising machines rapidly find low energy states [28]. In fact, higher-order oscillator Ising machines can satisfy more than  : Second-order versus higher-order networks when solving kSAT problems. a, The mean higher-order energy at the end of the simulation is plotted against the number of problem variables for hard instances of 3SAT problems. As the problem size increases, the difference in energy between the second-order oscillator Ising machines and the higher-order oscillator Ising machines increases. b, The mean percent of constraints satisfied at the end of simulation versus the problem size for 3SAT problems. c, The probability of satisfying all constraints for different problem sizes and models for 3SAT problems. d, The mean percent of constraints satisfied at the end of simulation versus the problem size for higherorder oscillator Ising machines for 3SAT problems. e, The mean time to satisfy 95% of constraints for higher-order Ising machines for 3SAT problems. d-e, Lines indicate different linear annealing schedules for the sub-harmonic injection locking coefficients, q i . In all plots, error bars represent the sample standard deviation computed over problem instances and trial simulations. f, Comparing resources and solutions of 5SAT and 7SAT problems to their 3SAT reductions. Reducing kSAT problems to 3SAT for k > 3 increases the number of variables and connections (left two columns). The 5th-order and 7thorder Ising machines find lower energy states corresponding to a greater fraction of constraints satisfied compared to the 3rd-order Ising machine (right two columns). 95% in less than one cycle for all problems.
Many studies on solving kSAT problems for k > 3, first use an efficient method for reducing the problem to 3SAT [41] and then focus on solving the resulting 3SAT problem. Here, we use a benchmark dataset of 5SAT and 7SAT problems [44] to assess this strategy for the higher-order oscillator Ising machine in terms of resource efficiency and solution quality (Methods 4.4). First, we find that the reduction to 3SAT increases the number of problem variables by one or two orders of magnitude, and there is approximately a 3 and 6 times increase in the number of clauses for 5SAT and 7SAT, respectively (left two columns in Fig. 3f ). Second, we observe that the direct solution of the 5SAT and 7SAT problems satisfy a greater fraction of constraints compared to solutions of corresponding 3SAT reductions (right column in Fig. 3f ). It would be interesting to compare the 5th-and 7th-order oscillator Ising machines to second-order oscillator Ising machines but we were unable to test second-order oscillator Ising machines on these problems due to the large number of auxiliary variables introduced via quadratization.

Discussion
Much of the existing literature on optimization with Ising machines have focused on second-order Ising networks. Such models were first proposed in [2] for solving constraint satisfaction problems. The authors in [2] originally proposed mapping a SAT problem to a higher-order polynomial but then applied quadratization to map to a second-order Ising model. Our first contribution is to directly compare the resource use of second-and higher-order Ising models for solving SAT problems. Defying common intuition, the comparison reveals that higher-order Ising machines are more resource-efficient than secondorder Ising machines for solving large combinatorial optimization problems. The resource efficiency of higher-order models results from the fact that no auxiliary variables are required and many combinatorial optimization problems map to polynomials which correspond to a very sparse higher-order interaction graph. Thus, the savings in higher-order models are in the number of Ising variables, as well as in the number of connections (Fig. 2).
Our second contribution is to build a resource-efficient higher-order Ising machine with coupled oscillators and test it on benchmark datasets of SAT problems. Motivated by other recent work [21], [36], [37], [45], we investigated the implementation of a higher-order oscillator Ising machine in a coupled oscillator network. Our model resembles the one in [45], but still differs in several ways. First, we use Hopf oscillators which include amplitude dynamics and capture the dynamics of oscillator hardware [46] more closely than the oscillators modeled by the Kuramoto model in [45]. Second, we introduce a form of annealing, specifically, the gradual increase of the sub-harmonic injection locking coefficient following a linear annealing schedule. Third, our model uses the simplest energy function resulting from the mapping method in [2] (Eqs. (1) or (2)), a sum of all constraint terms where each constraint term is a product of binary values. In principle, the mapping method specifies an entire family of valid energy functions in which the products in the constraint energy are raised by any positive exponent before summing them. For example, in [4], [45] the constraint terms are squared before summing. Our model choice results in gradient computations with the lowest possible complexity, and, moreover, achieves better solutions on the benchmark problems than a model with squared constraints (Methods 4.8).
Higher-order oscillator Ising machines converge to optimal or near-optimal solutions in very few cycles, and importantly, convergence time does not increase with problem size (Fig. 3e). In some practical cases, solutions are reached in less than one cycle. Further, higher-order oscillator Ising machines outperform second-order Ising machines in solution quality and in some cases find optimal solutions to Boolean constraint satisfaction problems. To our knowledge, this study is the first to report an Ising machine that finds optimal satisfiable solutions for the large 3SAT problems in the benchmark dataset (Fig. 3c).
It has to be emphasized that our study focuses on optimization methods with a basic Ising model whose only dynamic variables are the spin variables. These methods are extremely fast and resourceefficient, but they sometimes find only near-optimal solutions. Another type of Ising machine with higherorder interactions implements the Lagrange method [4], [47], [48], consisting of two types of dynamic variables, spin variables and Lagrange multipliers. In these models, each constraint term in the objective function is multiplied with a nonnegative variable, the Lagrange multiplier. If a constraint is unsatisfied, the corresponding Lagrange multiplier grows dynamically, until the constraint is satisfied [4], [47], [48]. In theory, the Lagrange models can find optimal solutions in polynomial time but the multipliers can grow exponentially large as a function of time [4]. Further, the time to solution in Lagrange models increases with problem size [4], [48]. The systematic comparison of Lagrange methods with higher-order versus second-order interactions is an interesting topic for future research.
The reported benefits of higher-order Ising machines, and higher-order oscillator Ising machines, in particular, are practically relevant because today many technologies exist for their realization. For example, higher-order interactions require the multiplication of the variables involved in the interaction. The multiplication of coupled electrical ring oscillator voltages can be implemented in the analog domain using existing CMOS technologies [49]. Further, the kth-order interactions of electrical oscillators can be implemented in log 2 (k) stages using a cascade of two-input multipliers or in one stage by a sequence consisting of element-wise log transform, summation, and anti-log transform. Another interesting technology is translinear electronic circuits which make use of the translinear principle [50]. Finally, existing methods for implementing real-valued analog higher-order interactions [51] may be modified for use in higher-order oscillator Ising machines.

Mapping optimization problems to higher-order Ising models
The energy function of the generalized Ising model, which includes higher-order interactions, is: Here, the real-valued variable J (k) represents the k-th order interaction between k spin variables and n is the total number of spin variables in the Ising model. The three groups of terms with 0th to 2nd order interactions on the RHS of (5) form the energy function of the traditional Ising model. Note that Eqs. (1) and (2) can be expanded and reduced into the form of (5).
In order to express the objective function of a combinatorial optimization problem as the energy function of an Ising model, binary variables in the optimization problem must be mapped to the spins of the Ising model. In this study, the transformation between problem variables, x i ∈ {0, 1}, and spins, s i ∈ {−1, 1}, uses the standard transformation: s i = 2x i − 1.

Equivalence of higher-order Ising energy formulations
It is easy to see that Eqs. (1) and (2) are equivalent. For constraint h, the corresponding setsC h or C h partition the state space. Therefore, any state, s, is an element of one of the two sets and we have: with Eq. (1) on the LHS and Eq. (2) on the RHS. The product terms evaluate to 1 when s = c and 0 otherwise. If s ∈C h , both sides equal 1, if s ∈ C h , both sides equal 0. Therefore, Eqs. (1) and (2) represent the same objective function and can be used interchangeably.

Derivatives of higher-order Ising energy functions
The partial derivatives of Eqs. (1) and (2) with respect to a complex variable, z i , can be efficiently computed as: and

Method for reducing k SAT to 3SAT
In this study, we also investigate a polynomial-time method [41] to reduce kSAT to 3SAT when k > 3. The method works as follows. Let ∨, ∧, and ∼ denote the logical OR, AND, and NOT operations, respectively. Consider a clause with 5 binary variables, ( . Introduce auxiliary variables y 1 and y 2 . Introduce new clauses and insert auxiliary variables as: The problem is 3SAT as no clause has greater than 3 variables. Reducing one 5SAT clause to 3SAT form results in 3 clauses and 7 variables. Consider a clause with 7 variables, ( . Introduce auxiliary variables y 1 , y 2 , and y 3 . Introduce new clauses and insert auxiliary variables: The last clause contains 4 variables so it has to be reduced further. Introduce auxiliary variables l 1 and l 2 . Introduce new clauses and insert auxiliary variables: The problem is 3SAT as no clause has greater than 3 variables. Reducing one 7SAT clause to 3SAT results in 6 clauses and 12 variables.

Excess resource use by different quadratization methods
Numerous quadratization methods have been proposed for reducing objectives with higher-order interactions to energy functions of second-order Ising energies [2], [17], [22], [24], [43]. In general, the number of auxiliary variables introduced by quadratization depends on the particular combinatorial optimization problem and the method of quadratization. In this study, quadratization was performed with the D-Wave Ocean software package 1 . With the quadratization method in D-Wave Ocean one can adjust the minimum energy gap, ∆E min , for controlling the tradeoff between excess resource use and computation performance of the resulting second-order Ising machine.
With the parameter choice of ∆E min = 1, one auxiliary variable is introduced per 3SAT clause, the same number as with some other quadratization methods [43]. Thus, the excess resource use of quadratization we report for this parameter choice generalizes to other methods in the literature. In addition, we also assess the resource use with parameter settings of ∆E min > 1 in D-Wave Ocean. These results are specific to the D-Wave Ocean quadratization method, but informative for exploring whether increased excess resource use could potentially close the performance gap between second-order and higher-order oscillator networks.

Benchmark datasets
We assess the performance of higher-order Ising machines on Boolean satisfiability (k-satisfiability, kSAT) problems, a well-known class of hard combinatorial optimization problems. Specifically, the 3SAT problems used in our experiments were obtained from the SATLIB collection [39] 2 . We selected instances of sizes 20, 50, 100, and 250 variables. The first sixteen instances were selected from each problem size to run the simulations. The dynamic variables in the oscillator networks were randomly initialized for each trial simulation. 64 trial simulations were performed for each instance.
To demonstrate the performance of higher-order Ising machines on 5SAT and 7SAT problems, we selected an instance of each problem from the 2018 SAT Competition [44] 3 . The 5SAT and 7SAT problems were also reduced to 3SAT using the method described in Section 4.4.

Oscillator model and simulation details
In higher-order oscillator Ising machines, each oscillator is represented by the complex Van der Pol or Hopf oscillator as described in Eq. (8): Here, ω i is the center frequency for the i th oscillator, λ i is a parameter determining the oscillator quality, and ρ i controls the degree of nonlinearity. In our simulations, the network coupling, r i (t), was the same for all oscillators and was held constant for the duration of the simulation. The center frequency was held constant at zero for all oscillators, ω i = 0 ∀i. The parameters λ i and ρ i were set to produce limit-cycle oscillations with unit amplitude. We used a linear annealing schedule, q i (t) = q max t t end . The phase quantization signal, l(z i ) is equivalent to sub-harmonic injection locking. We show this by representing each oscillator, z i , with a real and imaginary part (a i + ib i ). By adding the conjugate of z i to the dynamics, the real part grows and the imaginary part decays to zero. The solutions to the dynamics for each uncoupled oscillator including the limit-cycle dynamics, f (z i ), are a i = (h i + λ i )/ρ i The results reported in Fig. 3 were obtained using a parameter search to find the optimal values of λ, ρ, q max , and r. The best candidates were selected based on the lowest mean energy and the greatest mean probability of satisfying problem instances. The mean energy and percent of constraints satisfied were computed based on the final state of the network after simulation. The mean was computed across random network initializations for all trail simulations across problem instances within each problem size. The error bars in Fig. 3 represent the sample standard deviation. Integration of the dynamical system was performed using an adaptive step-size RK4/5 method. 4

Comparing higher-order constraint energy functions with different exponents
For a kSAT problem, the objective for clause h in our method (1) simplifies to: with c i = 1 if a literal is TRUE and c i = −1 if a literal is FALSE. Since E h (s) evaluates to either one or zero for all bipolar state vectors, s, an obvious generalization of the clause objective (9) is to exponentiate the RHS by a positive number. In [4], [45], the objective of kSAT problems with a higher-order energy function of this type was proposed, with the specific setting of the exponent set to a value of two: We compared the solution quality of higher-order oscillator Ising machines implementing objective (9) vs. (10). Our experiments included parameter optimization for each method, as described above (Section 4.7). Fig. 4 shows that networks based on (10) obtain worse solutions (with greater energies) and satisfy only a smaller percentage of constraints on benchmark 3SAT problems [39] compared to our method. The systematic analysis of exponent settings in the generalization of our method is left to future research. : Final energy on benchmarks 3SAT problems for the method proposed in this paper using Eq. (9) and the method using the square of the constraint energy (10) as in [4], [45]. a, Energy versus the number of variables for 3SAT problems. b, The percent of constraints satisfied versus the number of variables for 3SAT problems.