Combinatorial optimization with photonics-inspired clock models

NP-hard combinatorial optimization problems are in general hard problems that their computational complexity grows faster than polynomial scaling with the size of the problem. Thus, over the years there has been a great interest in developing unconventional methods and algorithms for solving such problems. Here, inspired by the nonlinear optical process of q-photon down-conversion, in which a photon is converted into q degenerate lower energy photons, we introduce a nonlinear dynamical model that builds on coupled single-variable phase oscillators and allows for efficiently approximating the ground state of the classical q-state planar Potts Hamiltonian. This reduces the exhaustive search in the large discrete solution space of a large class of combinatorial problems that are represented by the Potts Hamiltonian to solving a system of coupled dynamical equations. To reduce the problem of trapping into local minima, we introduce two different mechanisms by utilizing controlled chaotic dynamics and by dynamical formation of the cost function through adiabatic parameter tuning. The proposed algorithm is applied to graph-q-partitioning problems on several complex graphs. Physics-inspired algorithms are being developed to solve NP-hard problems while alleviating issues with scaling and trapping in local minima. Here, a method to search the global minimum of the Potts Hamiltonian with a photonic-inspired model is proposed.

C ombinatorial optimization problems such as graph partitioning and coloring problems appear in a wide range of applications 1 . These problems involve discrete search spaces and target finding the optimal solution from a finite set of objects. In computational complexity theory, it is known that many combinatorial problems are NP-hard 2 . Thus, solving large-scale combinatorial problems often means finding solutions with certain accuracy in a limited time rather than finding the globally optimal solution. On the other hand, it is well known that many important classes of combinatorial problems can be formulated as finding the ground state of lattice spin models including the Ising and Potts models 3,4 . In this regard, over the years, there has been a growing interest in simulating these models through quantum and classical physical systems for solving combinatorial optimization problems. These include superconducting qubits 5 , trapped ions 6 , and most recently photonic oscillator networks, including the coherent Ising machine 7,8 , the spatial Ising machine [9][10][11][12] , as well as polariton condensates and laser networks 13,14 .
Along different lines, numerical simulation of dynamical models governing physical systems has shown promise for solving combinatorial problems [15][16][17][18][19][20][21][22] . These schemes are built on direct simulation of nonlinear dynamical equations governing networks of interacting neurons with continuous degrees of freedoms that collectively evolve toward equilibrium points with binary states 17 . Despite their partial success, however, these methods suffer from two fundamental problems; (i) limitations of scaling the size of the network due to an increased computational burden, and (ii) trapping into local minima, which could be far from the global minimum of the governing cost function. Subsequently, there has been a continuous interest in developing novel physics-inspired algorithms that allow for finding good approximations of the solutions to large-scale combinatorial problems.
In this work, inspired by the dynamical equations governing a network of coupled q-photon down-conversion (q-PDC) oscillators, we introduce a new method for approximating the ground state of the planar Potts model. The proposed method is built on the dynamical evolution of a network of first-order nonuniform phase oscillators toward equilibrium states that generally represent the local minima of the Potts Hamiltonian. In order to mitigate the problem of trapping into the local minima, we utilize two different mechanisms that show significant improvement in minimizing the Potts Hamiltonian. The first approach is based on gradually forming the underlying cost function by tuning a parameter that controls the convexity at discrete phase values. The second approach is based on an adiabatic passage of the dynamical system through a chaotic regime which is enabled through numerical chaos that can be controlled by the timediscretization step size. Through examples, we show that the proposed minimalistic model in conjunction with a proper annealing mechanism enables an efficient solution of combinatorial optimization problems.

Results
The planar Potts model. The planar Potts model describes a network of interacting q-state spin degrees of freedom, where each spin can take an angle of θ = (2p − 1)π/q (p = 1, ⋯ , q) in the two-dimensional plane 4 . The interaction energy of a pair of spins (i, j) depends on their coupling strength J ij and the cosine of the relative angles (θ i − θ j ) according to H ij / ÀJ ij cosðθ i À θ j Þ. The planar Potts Hamiltonian is expressed as the sum of the pairwise interactions as 4 : which, is also referred to as the q-state clock model 23 . For q = 2, it reduces to the Ising Hamiltonian and for q → ∞ it approaches the XY Hamiltonian. Here, it is important to note the difference between the planar and the standard Potts models. The standard Potts Hamiltonian is defined as H ¼ À∑ <ij> J ij δðθ i ; θ j Þ, where, δ represents the Kronecker delta function. Clearly, the two models are identical for q = 2 and 3, while the planar Potts deviates from the standard model for q ≥ 4.
The Potts optimization problem. The Potts problem can be stated as finding the spin configuration that minimizes the Potts Hamiltonian. It can be readily viewed as coloring the vertices of a graph using q colors (q-coloring problem), where each vertex is represented by a spin and its color is represented by the spin state, and the graph edges weights are equal to the coupling strengths 3 . A number of other combinatorial optimization problems such as Max-q-Cut can be expressed in terms of the Potts problem. In this problem, we want to divide the vertices into q disjoint groups such that the total weight of edges between these groups is maximized. For q = 2, the Potts problem reduces to the Ising problem which represents the Max-Cut problem. These problems belong to the NP-hard computational complexity class meaning the solution space grows exponentially with problem size and it becomes prohibitively difficult to find optimum solutions to large-scale problems. For this reason, unconventional computing systems are being developed to find approximate solutions to large-scale problems in shorter times. It is worth mentioning that much of the research on unconventional methods and devices that target solving combinatorial optimization problems is focused on the Ising problem. These include the D-Wave's quantum annealing machine 5,24 , the coherent Ising machine 7,8 , Hitachi's CMOS Ising solver 25 , Fujitsu's digital annealer 26 , and Toshiba's simulated bifurcations machine 21 . In contrast, developing methods for solving the Potts problem offers an advantage in treating combinatorial problems that can be mapped more efficiently onto the Potts Hamiltonian.
The q-PDC oscillator. Recently, it was shown that the 3-state Potts model can be simulated with networks of 3PDC oscillators to optically solve certain optimization problems 27 . This is based on the tristable phase of the parametric spontaneous threephoton downconversion process 28 . The idea can be extended for building a classical nonlinear oscillator model for the q-photon downconversion (q-PDC) process, in which a photon of frequency qω is down-converted into q photons of frequencies ω. In this manner, we introduce a second-order nonlinear oscillator with complex amplitude a with the dynamics described as _ a ¼ Àγ l a þ gða Ã Þ qÀ1 À g s jaj 2ðqÀ1Þ a, where q = 2, 3, 4, ⋯ , and " * " shows complex conjugation. Here, g is the small signal oscillator gain, g s is the gain saturation coefficient, and γ l is the rate of intrinsic oscillator losses. The dynamics of this oscillator is governed by a Lyapunov function described as F 0 ¼ γ l jaj 2 À g q ða q þ ða Ã Þ q Þ þ g s q jaj 2q . Figure 1a schematically depicts the q-PDC oscillator and the associated virtual energy diagram for different q values. The corresponding single-oscillator Lyapunov functions are shown in Fig. 1b. These figures clearly indicate a phase transition as a result of increasing the gain g, which marks the onset of self-sustained oscillations with a finite amplitude. The equilibrium oscillator amplitude j aj obeys the nonlinear relation Àγ l j aj þ gj aj qÀ1 À g s j aj 2qÀ1 ¼ 0, which can be solved for various gain levels as depicted in Fig. 1c. According to these figures, the phase transition is second-order for q = 2 and first-order for q = 3, 4, ⋯ , where by increasing the gain a solution branch appears out of the blue sky 29 . For the purpose of optimization, a first-order phase transition causes additional complexities because it creates the undesired possibility of the death of oscillators 27 . However, as discussed later in this paper, by reducing the q-PDC oscillators to phase oscillators this problem will be effectively circumvented.
Next, we consider a network of linearly coupled q-PDC oscillators for q = 2, 3, 4, ⋯ . The equations governing the complex amplitudes of the oscillator network follows: where, a i represents the complex amplitude of the i'th oscillator, κ ij is the coupling coefficient between the i'th and j'th oscillators, and γ i = γ l + ∑ j |κ ij | is the damping factor of the i'th oscillator which is composed of the intrinsic oscillator loss γ l and external losses due to dissipative interaction with other oscillators ∑ j |κ ij | 27,30,31 . The relative strength of the self-oscillation g and coupling coefficients κ ij plays a critical role in the dynamics and accordingly in the success of the optimization process which will be discussed later in this paper.
The dynamical system described by Eq. (2) admits a Lyapunov function F , such that _ serves effectively as a cost function for the oscillator network. The Potts problem can be solved by numerically integrating the dynamical Eq. (2) which is equivalent to applying a gradientdescent minimization of the cost function of Eq. (3). Alternatively, one can start from the cost function F and apply any gradient-based technique for minimizing F over its continuous 2n-dimensional phase space 32 .
The nonuniform phase oscillators. The cost function of Eq. (3) embeds the planar Potts Hamiltonian, while a faithful representation of the Potts Hamiltonian requires uniform amplitudes for all oscillators 18,30,33,34 . To reduce the computational burden of numerical simulations, we consider uniform amplitudes for all oscillators, i.e., a i ! e iϕ i , which results in great reduction of the Fig. 1 From q-state oscillators to q-state Potts spins. a Schematic representation of a q-photon downconversion oscillator for q = 1, 2, 3, and 4. In this process, a pump photon of frequency qω is converted to q signal photons of frequencies ω (a p and a s represent the modal amplitudes of the pump and signal beams). b The associated energy landscapes above oscillation thresholds. c The phase transition associated with the onset of oscillations at the gain of g th . d The one-dimensional azimuthal potential governing the dimension-reduced phase oscillators. e The associated multi-state spin degrees of freedom.
governing 2n-dimensional Lyapunov function to an n-dimensional potential, F ! V, where V is a function of the phase degrees of freedom ϕ i , ⋯ , ϕ n (In transition from Eq. (3) to Eq. (4), a factor of − 2 is multiplied by the term involving g. This translates into a shift in the q-ary phase states which will not affect the results.): This relation is in essence the Lyapunov function governing a network of nonuniform phase oscillators 29 . This new potential carries the parameter g, which remains as a control parameter to adjust the strength of the self-oscillation of each oscillator compared to its interaction rate with the network. It should be noted that for q = 2 this relation reduces to the cost function introduced in ref. 35 and has been implemented with circuits for solving the Ising model 36,37 .
In the energy potential of relation (4), the first term operates as a penalty term that tends to enforce the evolution of each oscillator toward a q-ary state. The second term in isolation is equivalent to the XY Hamiltonian without the external field term. Therefore, assuming that the first term can successfully drive the oscillators to q-ary states, ϕ i = (2p − 1)π/q, p = 1, ⋯ , q, the potential reduces to the planar Potts Hamiltonian. However, this critically depends on the parameter g, which controls the strength of the penalty term.
Dynamical model and optimization. Considering the energy potential of Eq. (4), we can now write a dynamical model that guarantees descent motion toward the equilibrium points, ideally the global minimum. In doing so, we have the flexibility to deviate from the physics of the system and devise a dynamical model that results in faster and better numerical convergence to equilibrium points. Here, the straightforward approach is a gradient descent according to _ ϕ i ¼ À∂V=∂ϕ i , which results in: Alternatively, we can consider a pair of phase and angular velocity (ϕ i , ω i ) with dynamical equations _ ϕ i ¼ ω i and _ ω i ¼ À∂V=∂ϕ i À γω i , where γ is an effective friction coefficient. This choice is associated with the following second-order dynamics for the phase oscillators: To numerically find the equilibrium points, one can apply finite difference to either of Eq. (5) or Eq. (6) to iteratively update the phases at discrete time steps until reaching stable solutions. By applying different finite-difference formulas, different algorithms are obtained, while in general, momentum methods obtained from a discretization of Eq. (6) show faster convergence (please see "Methods").
In general, restricting the dynamics to a descent motion on the potential landscape of Eq. (4) could result in trapping into local minima. The problem of trapping into the local minima is a longstanding challenge in the optimization of non-convex functions. An interesting approach for avoiding this is to adiabatically increase the gain parameter 21,27,33 . This approach can be compared with quantum adiabatic optimization 38,39 . Although a fundamental difference is a lack of a classical counterpart for the quantum adiabatic theorem, which ensures a gradual evolution of the state of the system to the ground state of the final Hamiltonian. Here, the classical potential adiabatically evolves toward the desired combinatorial Hamiltonian, while the shallow local minima are less likely to arise during this process 27 . In the case of the phase oscillator networks discussed here, by gradually increasing the gain from zero, one can slowly guide the system into a valid (discrete) solution while staying inside the deeper regions of the coupling potential, possibly leading to a good approximation of the ground state solution. Clearly, it is important to increase the gain at an optimal speed that gives the network enough time to search for better solutions, while the caveat is that, slow tuning of the gain generally slows down the stabilization time. Figure 2a-c shows the landscape potential for an all-to-all connected network of n = 3 oscillators for three different gain levels. According to these figures, the potential deforms drastically by increasing the gain.
An alternative approach for bypassing the local minima is a controlled use of chaos. Chaos is by nature different from random Fig. 2 Adiabatic deformation and chaotic annealing mechanisms. For a 3-oscillator network (a), the potential landscape is shown in (b-d) for different gain levels g, and the dynamics is depicted in (e-i) for different discretization time steps Δ. These two parameters form the basis of the adiabatic deformation and chaotic annealing mechanisms that assist the optimization search process. a-c The local minima are shallower for lower than higher gains, which facilitates escaping the local minima when the gain is adiabatically increased. d-h The route to chaos as a result of increasing discretization time step, which creates the possibility of harnessing chaos for dynamical complexity in the optimization process. The chaotic map of a single oscillator for q = 2 (bi-stable oscillators) showing the oscillator's phase ϕ versus discretization time step is also depicted in (a).
fluctuations. However, the possibility of harnessing chaos for optimization has been suggested in previous works and surprisingly it is suggested that chaotic fluctuations can significantly outperform random noise 40,41 . A controlled transition of the oscillator network from chaos to fixed-point attractors allows for utilizing the exponentially fast searching capability of the network in its chaotic regime to largely avoid trapping into the local minima when reaching an equilibrium. It should be noted that an individual oscillator in Eq. (5), described by the dynamical equation _ ϕ ¼ g sinðϕÞ, is a first-order dynamical system and lacks the minimum dimensionality required for bifurcation and chaos. However, in principle, the dynamical model will be solved numerically, and as it is well known, numerical simulation of ordinary differential equations can induce chaos 42 . In this case, by using Euler's method, the phase of the oscillator at a time step t (k+1) can be obtained in terms of the previous step t (k) according to ϕ ðkþ1Þ It is straightforward to show that this simple iterative relation creates a chaotic map of the equilibrium phase ϕ versus the time step Δ. The chaotic map is depicted as an inset in Fig. 2 for a single oscillator with q = 2. As a result, by time-discretization, computational chaos can be incorporated in the system. Figure 2d-h depicts the dynamics of the three coupled oscillators for different values of the step size which clearly shows the chaotic dynamics for large step sizes. Tuning the parameter Δ allows a transition between equilibrium and chaos. Here, by starting with a large time step, the randomly initialized system evolves into a chaotic regime, and then by reducing the time step, the system makes an adiabatic passage from chaos to equilibrium.
Optimization results. We consider the Ising problem (q = 2) on an instance from the Biq-Mac library 43 , which involves a collection of quadratic binary programming problems. We chose the example of the "g05_100.6" graph, which is an unweighted graph with n = 100 nodes and an edge probability of 0.5 between any two nodes, containing a total of 2475 edges. Any other graph from this dataset essentially provides similar insights. Figure 3 shows the results of numerical simulation of the dynamical model of Eq. (5) with constant gain (a-c), adiabatic gain (d-f), and chaotic annealing (g-i).
Thefirstcolumndepictsthetime-evolutionoftheIsingHamiltonian (Eq. (1)) for 1000 random initial conditions. The distribution of the Fig. 3 Constant gain, adiabatic gain, and chaotic annealing methods for solving an Ising problem instance. The optimization results for the Ising problem on the "g05 100.6" graph with 100 nodes for constant gain (a-c), adiabatic gain (d-f), and chaotic annealing (g-i). The first column depicts the time-evolution of the Hamiltonian H ¼ ∑κ ij cosðϕ i À ϕ j Þ=2 (where κ ij is the coupling strength of the ith and jth oscillators and ϕ i is the ith oscillator's steady-state phase) for 1000 random initial conditions. The second column shows the distribution of the Ising energy H = ∑κ ij δ(ϕ i , ϕ j )/2 associated with the steady-state solutions. The third column exemplifies the dynamics of the phases ϕ for one simulation instant. For better visibility, the phase of each oscillator is shown with a different color. In (a-c), the gain is set to g = 2.5 and the step size is set to Δ = 0.01. In (d-f), the gain is linearly increased from g = 0 to g = 2.5 and the step size is set to Δ = 0.01. In (g-i), the gain is set to g = 2.5 and the step size is linearly decreased from Δ = 10 × 0.01 to Δ = 0.01. In the first column, the transient energy is calculated based on the XY Hamiltonian which transforms to the Ising Hamiltonian as the phases become binary after several iterations. These figures show the the minimum-to-maximum range of the energy distribution for all random initializations (light shaded area), the first standard deviation ('1 std dev') of the associated distribution (dark shaded area), and the average energy H ave (blue line) as a function of the iteration number.
Ising energy associated with the steady-state solutions is shown in the second column. The last column exemplifies the dynamics of the phases for one simulation instant. For adiabatic deformation of the energy landscape (d-f), the gain is linearly increased from zero. It is observed that at over 45% of the initial conditions lead to the bestknown solution (ground state), compared to 5% for the constant gain scenario (shown in panels (a-c)). By comparing the dynamics of the phases in panels (c) and (f), it is apparent that for the case of constant gain, the phases quickly settle into their final values. This is expected given the presence of many local minima that could involve the initial phase state in their basin of attraction. In this case, the dynamical system rapidly stabilizes in a local minimum. In the adiabatic gain regime, on the other hand, excursions continue to happen for longer simulation times, which indicates the evolving nature of the energy landscape.
For the case of the chaotic search (Fig. 3g-i), the gain is set to a constant value, similar to the first row, but the step size is reduced linearly, from Δt i = 10Δ 0 to Δt f = Δ 0 . According to Fig. 3h, 40% of initial conditions lead to the ground state. Here, the distribution of the steady-state energies turns out to be broader than the adiabatic gain case, and spanning over higher energy levels that are associated with local minima. Here, the rapidly evolving chaotic motion of the phases as depicted in Fig. 3i allows the system to explore distant parts of the phase space that otherwise would not be accessible.
Next, we apply the proposed method to solve the Ising problem on 100-node instances from all classes in the Biq-Mac library 43 .
This library involves nine classes of 100-node graphs, each containing ten instances. The results are shown in Fig. 4. In this figure, each panel shows the distribution of the steadystate energies for 1000 random initial conditions. The horizontal axis is the total number of steps, ranging from 10 steps up to 700. According to this figure, in almost all cases the top 20% of the ensemble reach the ground state with under 200 time steps.
In contrast with the Ising problem, to the best of our knowledge, the Potts problem (considering the case of q = 3) lacks a collection of benchmark problems with known solutions. Here, we generated sets of randomly connected 3-colorable graphs with different sizes and connection densities, with each set having 100 random instances. In general, average graph density or average degree of the graph is known to be a measure of the hardness for coloring problems, acting as an order parameter, with low-density graphs being generally harder to color 44 . Therefore, we focused on relatively low-density graphs with 500 and 1000 nodes and average degree of 9 (Fig. 5a, c) and 18 (Fig. 5b, d) (see "Methods" for details). In these simulations, the time-step size was set to a constant 0.01, while the gain is adiabatically tuned from 0 to 2.5. Each panel in Fig. 5 shows the distribution of the Potts energy associated with the steady-state phase of the oscillator networks versus the total number of time steps. The phases are initiated with uniform random values in the range − π ≤ ϕ i (0) ≤ π. Here, again, the top 20% of the solutions find the ground state in less than 200 steps.
all cases 50th percentile 20th percentile Fig. 4 Adiabatic gain method applied to the Ising problem. Solution to the Ising (q = 2) problem on 100-node graphs from all 9 classes in the Biq-Mac Library with known ground states. The average cut size C ave corresponding to the final average energy H ave is also shown on each panel. a-i The distribution of the equilibrium state energies H ¼ ∑κ ij cosðϕ i À ϕ j Þ=2 (where κ ij is the coupling strength of the ith and jth oscillators and ϕ i is the ith oscillator's steadystate phase) versus the number of iterations for 1000 random initial conditions. The horizontal axis shows the total number of steps (from 10 to 700). In all cases, the gain is linearly increased from g = 0 to g = 2.5, while the time-step size is Δ = 0.01.

Discussion
In summary, inspired by the q-photon downconversion parametric oscillator, we introduced a nonuniform phase oscillator that emulates a q-state spin in the two-dimensional plane. This provided an efficient route for searching for the ground state of the planar Potts Hamiltonian by numerically integrating the dynamical equations governing a network of phase oscillators that collectively evolve toward an equilibrium point. The dynamics of the oscillator network is restricted in the phase space to a Lyapunov function which is, in general, a non-convex function. Therefore, it is important to devise schemes that can prevent trapping into the local minima of the Lyapunov function.
Here, to reduce the problem of trapping into the local minima, we considered two different approaches based on the gradual formation of the cost function and by controlled passage through chaos. Our numerical results confirm that both schemes greatly improve the success of solving the optimization problem. The proposed method is ideally suited for parallel computing and application to large-scale optimization problems that can be formulated as the planar Potts model. This method is based on two parallel matrix-vector multiplications in conjunction with nonlinear function evaluations. In terms of the computational burden, the matrix-vector multiplication is the bottleneck of the proposed method. Therefore, speed-up can be achieved by parallel computing methods and systems.
It is worth stressing that for q ≥ 4, the clock Hamiltonian deviates from the standard Potts Hamiltonian. In this case, terms associated with unequal spin angles are penalized unevenly and depending on the relative angle of the two interacting spins. In this case, the clock model can still be used to approximate the ground state of the standard Potts model, while the approximation clearly deteriorates for larger q values.

Methods
For numerical calculations, a simple approach is to discretize the first-order dynamical model of Eq. (5) by using the first-order forward difference relation for the derivative term. This results in the Euler's method for solving ordinary differential equations which is based on a single-step updating of the dynamical variables in subsequent time steps k according to: where, Δ k = t (k+1) − t (k) is the time-step size, Φ ¼ ðϕ 1 ; Á Á Á ; ϕ n Þ t , and ∇V is the potential gradient, ∇V ¼ Àg sinðqΦÞ þ cosðΦÞ K sinðΦÞ À sinðΦÞ K cosðΦÞ; where, K is the connectivity matrix with elements κ ij and ' ⊙ ' shows the entryways vector product. In addition, we considered discretization of the second-order model of Eq. (6) by utilizing the central difference formula for the second derivative and the backward difference formula for the first derivative. This results in a two-step updating rule (known as the momentum method in numerical analysis): where, μ k = 1 − γΔ k and η k ¼ Δ 2 k . The convergence rate of the difference relations (9) can be controlled with the two parameters μ k and η k . In numerical simulations, for the adiabatic deformation approach, we used the two-step relation (9) by taking μ k = (k − 1)/(k + 2) and with a constant η k . For the chaotic annealing approach, we used the single-step relation (7), with a linearly decreasing η k .
The random graphs used for the 3-coloring problem were generated as described below. For a given number of nodes n, three partitions with n 1 , n 2 , and n 3 nodes were considered such that n 1 + n 2 + n 3 = n. Then, a full 3-partite graph was generated by connecting each node from one partition to all nodes of the other two partitions, while no edges exist in each partition. Next, the existing edges were removed at random with a probability of 0 < p 0 < 1, to create a random 3-colorable graph with the desired edge density. To obtain the energy distributions in each panel of Fig. 5, we first generated an ensemble of 10 2 random graphs with similar statistics, and in each case considered the a set of 10 2 random initial conditions. We then ran the simulations for all of the 10 4 combinations of random graphs and initial conditions, for each given total number of steps, and plotted the resulting energy distribution.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The code that supports the findings of this study is available from the corresponding author upon reasonable request. The edge weights κ ij are uniformly set to 1 for all the edges. For all panels, 100 random graphs are generated and each case is simulated with 100 random initial conditions to obtain the distribution of the energies H = ∑κ ij δ(ϕ i , ϕ j )/2 (where κ ij is the coupling strength of the ith and jth oscillators and ϕ i is the ith oscillator's steady-state phase) based on 10,000 data points for a given number of iterations. In all cases, the gain is linearly increased from g = 0 to g = 2.5, while the time-step size is Δ = 0.01.