Collective dynamics of phase-repulsive oscillators solves graph coloring problem
Aladin Crnkić, Janez Povh, Vladimir Jaćimović, Zoran Levnajić
aa r X i v : . [ n li n . C D ] M a r Collective dynamics of phase-repulsive oscillators solves graphcoloring problem
Aladin Crnki´c , Janez Povh , Vladimir Ja´cimovi´c , and Zoran Levnaji´c ∗ Faculty of Technical Engineering, University of Biha´c, Ljubijanki´ceva, bb., 77000 Biha´c, Bosnia andHerzegovina Faculty of Mechanical Engineering, University of Ljubljana, Aˇskerˇceva cesta 6, 1000 Ljubljana,Slovenia Faculty of Natural Sciences and Mathematics, University of Montenegro, Cetinjski put, bb., 81000Podgorica, Montenegro Complex systems and Data science Lab, Faculty of information studies in Novo mesto, Ljubljanskacesta 31A, 8000 Novo mesto, Slovenia Department of Knowledge Technologies, Joˇzef Stefan Institute, Jamova cesta 39, 1000 Ljubljana,Slovenia
Abstract
We show how to couple phase-oscillators on a graph so that collective dynamics ‘searches’for the coloring of that graph as it relaxes towards the dynamical equilibrium. This trans-lates a combinatorial optimization problem (graph coloring) into a functional optimizationproblem (finding and evaluating the global minimum of dynamical non-equilibrium potential,done by the natural system’s evolution). Using a sample of graphs we show that our methodcan serve as a viable alternative to the traditional combinatorial algorithms. Moreover, weshow that with the same computational cost our method efficiently solves the harder problemof improper coloring of weighed graphs.
Despite the explosion of modern computing power, problems of combinatorial opti-mizations remain a formidable algorithmic challenge. Drawing upon the theory ofcollective dynamics of phase-oscillators on graphs (networks), we show how to con-vert a combinatorial optimization problem (graph coloring) into a functional min-imization problem (finding and evaluating the global minimum of a non-negativefunction on a bounded domain). While this does not (necessarily) offer a sim-pler solution to the graph coloring problem, it lend it to the whole new world ofwell-developed methodologies – those of functional optimization.
Problems of combinatorial optimization pervade all walks of science [1, 2]. This is particularlytrue in modern data-driven era, as analyzing the data becomes much harder than obtaining thedata. Combinatorial optimization problems revolve around finding the optimal value of somecriteria function on a discrete set. That set is usually very large, typical example being theset of all combinations or variations of some elements in a given order. Classical combinatorialoptimization problems include linear and quadratic assignments, bin packing problem, vehiclerouting problem, graph partitioning, etc. [3]. ∗ To whom the correspondence should be addressed: [email protected]
The graph coloring problem
We begin be precisely formulating the (vertex) graph coloring problem as the combinatorialoptimization problem. Let G = ( V, E ) be a non-directed graph with vertex set V and edge set E . A K -coloring of vertices V is a function ϕ : V → { , , . . . , K } , which maps adjacent verticesinto different numbers, i.e., ϕ ( u ) = ϕ ( v ) for all uv ∈ E . Graph is called K -colorable if and onlyif there exists a K -coloring of it vertices. The minimal K , such that there exists a K -coloring iscalled chromatic number of graph G and is denoted by χ ( G ) [37]: χ ( G ) := min { K | ∃ K -coloring of vertices of G } . (1)Finding the chromatic number of given graph is referred to as graph coloring problem , and is aclassical problem of combinatorial optimization. There is a vast amount of theoretical results forthis problem on different families of graphs and on the problem’s extensions [38]. Graph coloringserves as a theoretical model for several practical problems, including timetable scheduling [39]or frequency assignment problem [40].The problem can be formulated in the context of optimization theory in several differentways [41, 42, 43]. Here we present the most basic version, so-called assignment-based ILPmodel [43]: χ ( G ) = min P Hi =1 w i s.t. P Hi =1 x vi = 1 , ∀ v ∈ Vx ui + x vi ≤ w i , ∀ uv ∈ E, i = 1 , . . . , Hx vi , w i ∈ { , } , ∀ v ∈ V, i = 1 , . . . , H
In this model H is a fixed upper bound for χ ( G ) obtained by applying some of the manytheoretical upper bounds or heuristic algorithms. Binary variables x vi represent the assignmentof vertices to colors. Each vertex must be assigned exactly to one color (the first constraint)while two adjacent vertices can not have the same color, which is controlled by binary w i (thesecond constraint).Graph coloring problem is simple for some special families of graphs. Trivial examples are K -partite graphs which are K -colorable. For the complete graph on N vertices K N is trivially χ ( K N ) = N . For perfect graphs, the chromatic number equals to the size of the largest cliqueand this number can be efficiently computed by semi-definite programming. This follows fromthe definition of perfect graphs and by the famous Sandwich theorem [44], while the Strongperfect graph theorem simply characterises these graphs [3].As for solving the graph coloring problem, exact algorithms that solve it to optimality startwith one of the mathematical programming formulations and try to feed it into appropriate non-linear programming solver, if the problem is small enough and if the there is a solver available forthat formulation [45]. For larger graph sizes, linear and semi-definite programming relaxationsare used within Branch and Bound framework. However, for practical needs, a vast amount ofheuristic algorithms is available [46, 47]. Coming back to phase-oscillators, we begin explaining the contribution of this paper via il-lustrative toy-example. Consider a graph with N nodes defined by the adjacency matrix A .Phase-oscillators specified by phases ϕ i are assigned to its nodes, each with degree k i . Oscil-lators’ frequencies are identical so we set them all to ω i = 0 without loss of generality. Thecoupling strength is ε = −
1. The equation for this phase-repulsive graph of oscillators reads:˙ ϕ i = − (1 /k i ) X j A ji sin( ϕ j − ϕ i ) . (2)3f our graph had only N = 2 oscillators (nodes) connected by a link, dynamics would push theirphases away from each other to the opposite values. Eventually, the phase difference betweenthem would reach π , since this is the stable equilibrium for this system [29]. Assume nowthat our graph is . Dynamics pulls the phases of the connected oscillators apart, i.e.,‘stretches’ the phase differences along each link, trying to reach the maximum stretch value π .Since the graph is 2-partite, this equilibrium is easily attained. Eventually, oscillators in twogroups will have the opposite phase values and the graph will be perfectly anti-synchronized.An illustration of this situation is shown in Fig. 1a.Figure 1: (a) Illustration of a simple 2-partite graph (network): link connect only nodes in onegroup with the nodes in the other group. Groups are shown in grey. (b) Illustration of a simplegraph that is not 2-partite, but it is 3-partite.But now, inverting this problem, Eq. 2 can be used to check whether a given graph is 2-partite. Once the dynamics of Eq. 2 reaches equilibrium, one looks at the phase differences alonglinks: if they are all π the graph is 2-partite. Actually, consider a graph that is not 2-partite.For simplicity, let it have 3 nodes connected in a triangle, like in Fig. 1b. It is a 3-partite graph,since nodes can be grouped in three groups with links going only between the groups (but not intwo groups). The dynamics of Eq. 2 will no longer be able to ’stretch’ the phase differences alonglinks to π . Instead, at least some links will be “frustrated” or “squeezed” to a phase differenceless than π . This will happen for any choice of initial phases, indicating that the underlyinggraph is not 2-partite.To articulate this dynamical situation we introduce frustration f ij = f ij ( ϕ j − ϕ i ) along thelink i − j ( A ij = 1) as: f ij ( ϕ j − ϕ i ) = 1 + cos ( ϕ j − ϕ i ) . (3)Frustration captures how “squeezed” is a link, f ij = 0 means that link i − j is not frustrated(phase difference is π ). If we picture the oscillators’ interactions as springs, frustration can bethought of as elastic potential energy due to spring being squeezed [29]. In fact, gradient of f ij yields the sinusoidal interaction in Eq. 2, so f plays the role of dynamical non-equilibriumpotential . We define the total frustration F of a graph as the sum of frustrations along theindividual links: F = X i,j A ij f ij ( ϕ j − ϕ i ) , (4)so that F = 0 guarantees that f ij is zero along all links. The graph in Fig. 1a will always have F = 0, whereas the one in Fig. 1b will always have at least some frustrated links and hence F >
0. This shows that former graph is 2-partite, while the latter is not. Eq. 2 can be used toconfirm 2-partitness, but not to confirm 3-partitness.So, can we develop a dynamical model inspired by Eq. 2 that could confirm K -colorabilityof a graph for any K (not just for K = 2)? Such model would offer a dynamical solution to Frustration can be equivalently defined for phase-attractive case, but it is always zero since for identicaloscillators full synchronization is the only equilibrium state. In contrast, for the phase-repulsive case graphtopology is intimately related to frustration [29] dynamically , just bynaturally relaxing to its equilibrium. Note the analogy with searching for the ground state inquantum systems [16] (quantum algorithms here serve only as inspiration and have no relevanceotherwise).Equipped with above preliminaries, we now present exactly such a model. We are givena non-weighted and non-directed graph and set out to find a dynamical model that checks its K -colorability for any K . We need the specific coupling functions and frustration function foreach K . For K = 2 we have Eq. 2 and Eq. 3: a graph is 2-colorable if for at least one choice ofinitial phases we obtain F = 0.We define the K -frustration f Kij along the link i − j for K ≥ f Kij ( ϕ j − ϕ i ) = 1+ c K h ( K −
1) cos( ϕ j − ϕ i )+( K −
2) cos 2( ϕ j − ϕ i )+ · · · +cos( K − ϕ j − ϕ i ) i . (5)Constants c K are fixed to have f Kij ( πK ) = 0 for every K . Each f K is a trigonometric polynomialof the degree K −
1. For K = 2 we have f ij = 1 + c cos( ϕ j − ϕ i ). Requesting f ij ( π ) = 0 we get c = 1 and recover Eq. 3.We define the total K -frustration F K similarly to Eq. 4 as the sum of K -frustrations alongthe individual links F K = X i,j A ij f Kij ( ϕ j − ϕ i ) . (6) F K is the non-equilibrium potential for the dynamical model behind Eq. 5. We have:˙ ϕ i = − k i ∂F K ∂ϕ i = − k i X j A ij ∂∂ϕ i f Kij ( ϕ j − ϕ i ) = − k i X j A ij C K ( ϕ j − ϕ i ) . (7)For coupling functions C K we get: C K ( ϕ j − ϕ i ) = c K h ( K −
1) sin( ϕ j − ϕ i )+2( K −
2) sin 2( ϕ j − ϕ i )+ · · · +( K −
1) sin( K − ϕ j − ϕ i ) i , where the constants c K are fixed as before. Since c = 1, Eq. 7 for K = 2 reduces to thephase-repulsive graph of Kuramoto oscillators Eq. 2.For higher values of K , in addition to the first harmonic like in Eq. 2, oscillators are alsocoupled via higher harmonics. Consider the example of K = 3. We have: f ij ( ϕ j − ϕ i ) = 1 + 43 cos( ϕ j − ϕ i ) + 23 cos 2( ϕ j − ϕ i ) , ˙ ϕ i = − / k i X j A ji h sin( ϕ j − ϕ i ) + sin 2( ϕ j − ϕ i ) i . (8)For phase differences ϕ j − ϕ i = π or ϕ j − ϕ i = π we get f ij = 0 and hence ˙ ϕ i = 0. Incontrast, for all other phase differences we have ˙ ϕ i = 0. Hence, above dynamical model looksfor equilibrium states where the phase differences along links are integer multiples of π . Backto the graph from Fig. 1b, under dynamics governed by F its three nodes will ultimately havethree different phase values, equidistant on the circle, and separated by π . In this situation wehave F = 0, indicating that this graph is 3-colorable. Meanwhile, F = F will stay positive,since graph is not 2-colorable (every K -colorable graph is also ( K + 1)-colorable, but not viceversa).In general, any graph is K -colorable if and only if the global minimum of F K is zero and thesmallest such K ≥ χ ( G ). In fact, if for some arrangement ofequilibrium phase values around the nodes we find F K = 0, then there are exactly K different5quilibrium phase values. This means some nodes will have common phase values, but neverthose that are connected. Instead, phases of the connected nodes will always be integer multiplesof πK apart (if they were not, this would not be an equilibrium state, and F K would be positive).Identifying K colors as K equilibrium phase values, we have that such graph is K -colorable.This is the core result of our paper and it is rigorously proven in the Appendix.To confirm that a graph is K -colorable one needs to run the dynamical model Eq. 7 frommany initial phases and examine the final equilibrium states, searching for the one with F K = 0.Simplest way to do this is to track F K ( t ) (in addition to several ‘tricks’ that we discuss later).Note that F K always has a global minimum: the question is, if it is zero. Thus, checking the K -colorability of a graph is translated into evaluating the global minimum of a real function F K ,which is done via natural evolution of the collective dynamics Eq. 7. We have thus translateda combinatorial into a continuous optimization problem: rather than permuting arrangementsof node colorings (combinatorial problem), we run the Eq. 7 by integrating ordinary differentialequations (continuous problem).However, if F K = 0 is not found, this still does not prove that graph is not K -colorable. F K is a trigonometric polynomial whose functional properties depend on (and in fact, capture) thetopology of the underlying graph. Therefore, in addition to the global minimum, such functionwill in general have many local minima where F K >
0. Dynamics could easily get “stuck” inone of them, and confuse it for the global minimum, leading to (potentially false) conclusionthat graph is not K -colorable. Dynamics can miss the true global minimum for several reasons,e.g. because its basin of attraction is too small. This is a notorious problem in functionaloptimization.Furthermore, there are situations where the minimization converges not to a point, but toa critical set of bigger dimensionality, for example a limit cycle. However, this does not changethe principle of our method: we look for F K = 0, regardless of the dimensionality of the set onwhich this occurs. In fact, larger the set with F K = 0 is, easier is for convergence process to findit. On the other hand, convergence can become stuck more easily in a local minimum of biggerdimensionality (than in a point). Partial remedy for this is in ‘tricks’ that we discuss later. An important remark.
Once the non-negative function F K is defined, its local and globalminima are what they are (zero or otherwise), regardless of whether our dynamics Eq. 7 findsthem or not. So, running the dynamics Eq. 7 is not the only way to evaluate the global minimumof F K . There are numerous approaches in functional optimization for evaluating the minima ofreal functions, such as gradient descent, differential evolution, etc [48]. In fact, some of themcould be more efficient than the natural evolution of Eq. 7. In other words, once the combina-torial problem of K -colorability is converted into the optimization (evaluation) problem, it canbe approached either via Eq. 7 or via any of the standard methods in functional optimization. To illustrate our result, we construct a simple graph with 8 nodes and chromatic number 4. Werun the dynamics of Eq. 7 with coupling function C , whose corresponding frustration is F .Initial phases are set at random, but instead of picking them from uniform distribution on thecircle, we use von Mises distribution [49] (it can roughly be thought of as a Gaussian on thecircle, where κ ≥ κ = 0). The effect of this is that initially all phases are close to each other, so the system isvery frustrated, which gives it a bit of “push” and makes phases diverge more rapidly.In Fig. 2 we show two realizations of dynamical evolution on the graph via four snapshots(last snapshot for t = 15 is the equilibrium state). In top panel the dynamics reaches the statewith F = 0, which confirms that graph is 4-colorable. There are indeed only 4 different final6hase values and no connected pair of nodes is colored the same. In bottom panel the dynamicssettles into a state with F >
0, since some links remain frustrated and there are more than 4different equilibrium phase values. This is a local minimum, in the sense that small perturbationswill not “kick” it out of it. This shows how not all runs finish in the global minimum, regardless t = = = = phases π π π frustration Figure 2: Two realizations of evolution for illustrative graph with 8 nodes, shown in four timesnapshots as indicated above. Both evolutions start from initial phases randomly chosen fromvon Mises distribution on the circle with κ = 10 [49]. The phases of the nodes are shown bycolors and frustration of links in gray-scale (see two colorbars on the right). Top panel is anexample of dynamics that reaches F = 0, while bottom panel is an example of dynamics thatsettles in a local minimum with F > F = 0 was much narrower? In thatcase it could happen that no runs finish in the global minimum. This would falsely lead us tobelieve that graph is not 4-colorable. Establishing whether some minimum is local or global isa notorious problem in functional optimization. To investigate its repercussions here, we trace F K ( t ) starting from 50 initial phases for above graph and report the results in Fig. 3. Situationin earlier Fig. 2 corresponds to Fig. 3b ( K = 4). Only a fraction of runs finish in F = 0 (redcurves), and the rest finishes in F > F = 0, assuming there areany. To put this in context, we re-do the same runs with K = 3 are show the plots in Fig. 3a. F a ) t F b ) t F c ) t Figure 3: (a) Plots of F ( t ) for 50 runs of Eq. 7 with C for the same graph as in Fig. 2. (b)The same for C and F , two examples on this dynamics are shown in Fig. 2. (c) The same for C and F . Plots are shown until t = 15 for consistency with Fig. 2. Red curves are runs thatfinish in F K = 0 and blue curve are runs that finish in F K > F = 0, since graph is not 3-colorable. However, the same7ituation could happen for K = 4 unless we run enough realizations. In contrast, in Fig. 3c weshow the same for K = 5: almost all runs lead to F = 0, since graph is easily colorable with 5colors. Confirming the colorability for higher K is easier, since F K touches zero in more points,despite F K being a higher-order trigonometric polynomial for higher K . As we get closer to theminimal K for which the graph is K -colorable, less and less initial phases lead to F K = 0, sincethe basin of attraction for states with F K = 0 shrinks. That is why pinpointing the chromaticnumber of a graph is challenging, as illustrated in Fig. 3.We ran many more simulations for graphs of varying size and chromatic numbers. Overall,we found that our simulations correctly identify the chromatic numbers obtained via traditionalcombinatorial algorithms. However, as graph size grows beyond N = 50, the number of necessaryruns increases rapidly. In fact, the dimensionality of space on which F K is defined is N . As F K becomes higher and higher dimensional, pinpointing the basin of attraction for F K = 0 becomesharder and harder. But still, if the chromatic number is (relatively) low, the corresponding F K is a low-order trigonometric polynomial, so our dynamical graph coloring works.On the other hand, except taking the initial phases from specific distributions, there are otherenhancements (or ‘tricks’) one could use. One of them is to introduce kicking in the oscillatordynamics, i.e., make Eq. 7 a stochastic differential equation by adding a noise term on its RHS.The effect of this is that evolution can get kicked out of a local minimum, possibly finding theglobal minimum later. This is similar to the effect that temperature has in simulated annealing.We confirmed that this enhancement indeed improves the overall performance of our method,but it still does not solve the challenge of large graphs and high chromatic numbers. Our approach is actually useful for a harder problem known as improper graph coloring [50, 51].Here, one looks into colorings with number of colors smaller than graph’s chromatic number.Therefore, one allows for some connected nodes to be colored the same, but looks for thearrangement of colors that minimizes the number of such pairs. In this improper sense, anygraph is colorable with any K number of colors, but the question is, how to minimize thenumber of connected nodes having the same color. If this number is zero for some K , the graphis K - colorable in the proper sense, and we have χ ( G ) ≤ K . Several variations of this problemhave been studied in the literature, although not very frequently, with diverse motivations, forexample related to scheduling and timetabling problems [52].To make our analysis even more general, in this section we face the problem of impropercoloring of weighted graphs. In a weighted graph, links (edges) are not all of equal weight(“thickness”), but can be stronger or weaker, modeling the strength of interactions betweenthe nodes. Consider a non-directed weighted graph with N nodes, described by the symmetricweighted adjacency matrix W ij . The entry W ij in this matrix denotes the (non-negative) weightof the link i − j . Given an integer K ≤ N , the problem is to color the nodes of this graph in away to minimize the sum of weights of links that connect nodes that are colored the same. Inother words, we want to minimize the sum of W ij such that nodes i and j are colored the same(if i and j are not connected then of course W ij = 0). Traditionally, improper graph coloring istreated via similar combinatorial algorithms, which in weighted case have an additional layer ofcomplexity.As before, we assign phase-oscillators of identical frequencies ω = 0 to graph nodes andconsider a tentative coloring number K . K -frustration f Kij along i − j remains defined by Eq. 5,but the total K -frustration F KW is now a weighted sum of individual link frustrations: F KW = X i,j W ij f Kij ( ϕ j − ϕ i ) . (9)8he dynamical equations are obtained the same way:˙ ϕ i = − k i ∂F KW ∂ϕ i = − k i X j W ji ∂∂ϕ i f Kij ( ϕ j − ϕ i ) = − k i X j W ji C K ( ϕ j − ϕ i ) , (10)where the coupling functions C K for each K are the same as before. Our main result fullygeneralizes: if the system reaches the global minimum of F KW , then in that state there areexactly K different equilibrium phase values on the nodes, separated by integer multiples of πK .The key difference is that now even in the global minimum we have F KW >
0, since some linkswill be frustrated, as their nodes will have the same phase value (i.e. be colored the same).However, the sum of weights W ij of such links will be minimal, as we prove in the Appendix.Hence, the minimization of F KW yields the improper coloring of a weighted graph using exactly K colors. As in the non-weighted case, we can use either the dynamics itself or a range of otherfunctional minimization methods.There are several differences with the non-weighted case. First, a graph is now K -colorablefor any K , the only question is, how frustrated it will be. Our results says that the dynamicsevolves to ( i ) minimize the total frustration F KW , and ( ii ) squeeze the remaining frustrationpreferentially into weaker links. We have again translated a combinatorial into a continuousoptimization problem, it is just that now we cannot recognize the global minimum via F KW = 0,but via K different final phase values. Second, the fact that coloring is improper does notspecify how many connected node pairs are colored the same (at least one), but only that thejoint frustration of their links is minimal. So, we are not minimizing the number of frustratedlinks, but the total frustration F KW itself. Third, in weighted case dynamics can also get stuck inthe local minima. For any K there is a guaranteed improper coloring corresponding to minimaltotal frustration, but one might again resort to “tricks” such as special choice of initial phasesor stochastic evolution to find that minimum more efficiently. On the other hand, local minimacan now be interpreted as Nash equilibria in the system: let e Ki be the K -frustration of thenode i defined as the sum of K -frustrations of all its adjacent links e Ki = P j W ij f Kij . Collectivedynamics can now be seen as each node i seeking to minimize its e Ki . Pairs of nodes connectedby weak links will be more successful in this.Next we show the performance of our dynamical method for improper coloring. For sim-plicity we use a fully connected graph (clique) with N = 6 nodes (this is a non-trivial problem,a variation of the channel assignment problem for wireless graphs known from telecommunica-tions). We randomly assign weights to its 15 links and run the dynamics of Eq. 10 with couplingfunctions C , C , C , C , computing the corresponding frustrations F W , F W , F W , F W . Theresults are shown in Fig. 4, three snapshots for each example. For C , as expected, the dynamicssettles into a state with zero total frustration and with graph (properly) colored with K = 6colors. For C , as the graph tries to get colored with 5 colors, at least one pair of nodes mustbe colored the same. The dynamics reaches a state with just one link frustrated, the weakestone. For C , the dynamics settles into a state where two weakest links are frustrated. For C ,graphs is eventually colored with 3 colors. Four links remain frustrated and we verified thattheir joint weight is the smallest possible. So, as expected, in all cases our dynamical coloringreached the equilibrium with minimal total frustration squeezed into weakest links. In caseof a complex topology (not clique), the minimization into weakest links would interplay withtopology, making it more difficult to see that algorithm works as expected.To examine the impact of the local minima, we repeat in Fig. 5 the analysis from Fig. 3.We track F KW ( t ) between t = 0 and t = 15 for 50 choices of initial phases for above graph for K = 6 , , ,
3, in correspondence to four panels in Fig. 4. For K = 6 most initial phases lead to F W = 0, but not all: since the graph is weighted, despite having 6 colors available, the graphsometimes fails to “stretch” to zero total frustration. In such equilibria, the graph is not colored9 = F t = =
15 t = F t = = phases π π π F F frustration Figure 4: Four realizations of evolution for weighted fully connected graph with N = 6 nodes,shown in three snapshots each, as indicated. Link weights are chosen randomly between 1 to 5.Initial phases are randomly selected from von Mises distribution with κ = 2. The phases of thenodes are shown by colors and frustration of links in gray-scale, as in Fig.2. One realization ofevolution for each coupling function C , C , C , C and the corresponding frustration F W , F W , F W , F W is shown in each panel, as indicated. F a ) t F b ) t F c ) t F d ) t Figure 5: (a) 50 realizations of F W ( t ) for 50 runs of Eq. 10 with C for the weighted graphfrom Fig. 4. (b) The same for C and F W . (c) The same for C and F W . (d) The same for C and F W (all in correspondence with four panels in Fig. 4). Plots are shown until t = 15 forconsistency. Red curves are runs that finish in zero total frustration and blue curves otherwise.by 6 colors. For K = 5 , K colors decreases with K . For each of these values of K there are many local minima,so the global minimum is somewhat less distinguishable, and has to be identified by actuallychecking if the graph is colored (although improperly) with K distinct colors. We run additionalsimulations for larger weighted fully connected graphs and found the same scaling propertiesas in the non-weighted case. Global minima are increasingly harder to find as the graph sizeincreases. Enhancements such as kicking (stochastic evolution) have the same effect as before –they improve the search, but does not solve the problem of local minima. Using the analogy with repulsively coupled phase-oscillators we found a way to convert a com-binatorial optimization problem (proper and improper graph coloring of non-weighted andweighted graphs) into the optimization problem of evaluating the global minimum of a realfunction. This function is the non-equilibrium potential for the collective dynamics of oscillatorscoupled through the edges of the graph (network). This dynamics can be used as a “natural” wayto search for the minimum, or alternatively, one can resort to the rich ensemble of methods fromfunctional optimization. This conversion does not offer an immediate solution to the problem,10ut it allows us to approach the problem using completely different methodology. Examiningin what situations this might lead to a more efficient or precise solution of the original problemremains the matter of future work, potentially extending to other discrete problems.There are two ways of comparing our result to the state of the art. First, one could compareour collective dynamics as a method of evaluating the global minima of real functions to theexisting optimization methods. We did not find this of interest, since this means comparingfunctional optimization methods among them, and for that there already exists abundant liter-ature. Second and more interesting way is to examine whether solving graph coloring problemthis way would be better than solving it via standard combinatorial algorithms. We did lookinto this, and as already stated, for all examined graphs we found the same chromatic numbers.In no case we found the chromatic number bigger or smaller than indicated by traditional al-gorithms. We checked this for 30 graphs with sized up to N = 100. Another independent wayof making this comparison is to look at the computational cost involved in two cases. However,this comparison strongly depends on the computing architecture, since the algorithms are ofquite different nature (permutations of color arrangements vs. integrating a differential equa-tion using e.g. Runge-Kutta integrator). An interesting direction of future work might be tofind a computationally cheaper integrator for Eq. 7 [53].A different issue revolves around scaling up the problem: how about coloring a graph ofsize N = 1000? This represents a challenge for both combinatorial and our method. While theformer suffer from combinatorial explosion, our method comes down to finding a minimum of1000-dimensional real function, which is far from trivial. However, while cases of small graphscan in practice be easily treated by any approach, the true challenge lies in coloring large graphs,which is where our approach might make a difference. This remains the core direction of futurework. Still, the degree of trigonometric polynomial in Eq. 5 is K −
1. Hence, checking thecolorability for small K is relatively easy, since low-level trigonometric polynomial is not very“wavy” and the minima is easier to find, even for large graphs. This might be a competitiveadvantage of our methods, since checking for these minima could beat the existing combinatorialoptimization approaches. In contrast, when checking the colorability for high values of K , thedegree of trigonometric polynomial is higher, which makes the function more wavy and hencethe basins of attraction for local and global minima are steeper and harder to find. On the otherhand, given the purely trigonometric nature of any F K , it would be interesting to do frequencyanalysis on it.We close the paper with another potentially interesting idea. Consider a graph G withunknown chromatic number χ ( G ) = K and some very large tentative K , much larger than K (say K = N ). F K is zero on a potentially very large domain within R N , since the graph iseasily K -colorable for this K ≫ K . We call this domain support of F K . Consider now K − F K − is a large domain, although presumably smaller than support of F K . Continuing on, one expects that support of F K − is even smaller, but larger than thesupport of F K − , and so on. Eventually, the support of F K is the smallest such non-emptydomain, since the support of F K − is empty. What is the geometric relationship between thesedomains? For example, could they be nested subsets of each other? If so, one could be able toextrapolate K by examining the geometric process of how these domains (supports) shrink as K deceases. The value of K just before this domain becomes empty is the chromatic number K , obtained without optimization. Acknowledgements.
Work supported by the Slovenian Research Agency (ARRS) via pro-grams P1-0383 and P2-0256, projects J5-8236, J1-8155, N1-0057 and N1-0071, and by EU viaMarie Sklodowska-Curie Grant project 642563 (COSMOS). Part of it has been completed duringthe STSM of the third author (VJ) supported by the COST action CA15140. Thanks to col-leagues Alexander Yurievich Gornov, Peter Koroˇsec, Ljupˇco Todorovski, and Riste ˇSkrekovski11or very useful suggestions and feedback.
References [1] L. A. Wolsey and G. L. Nemhauser,
Integer and combinatorial optimization (John Wiley& Sons, 2014).[2] C. H. Papadimitriou and Kenneth Steiglitz,
Combinatorial Optimization: Algorithms andComplexity (Courier Corporation, 2013).[3] A. Schrijver,
Combinatorial optimization: polyhedra and efficiency (Springer Science &Business Media, Vol. 24, 2003).[4] M. Baghel, S. Agrawal, S. Silakari. ”Survey of metaheuristic algorithms for combinatorialoptimization.” International Journal of Computer Applications 58, 19 (2012).[5] B. Korte and J. Vygen,
Combinatorial Optimization: Theory and Algorithms (SpringerScience & Business Media, 2007).[6] G. J. Woeginger, ”Exact algorithms for NP-hard problems: A survey,” in
CombinatorialOptimization — Eureka, You Shrink! , edited by M. J¨unger, et al. (Springer, Berlin Heidel-berg, 2003), pp. 185-207.[7] S. Salcedo-Sanz, ”Modern meta-heuristics based on nonlinear physics processes: A reviewof models and design procedures,” Phys. Rep. 655, 1-70 (2016).[8] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, ”A quantumadiabatic evolution algorithm applied to random instances of an NP-complete problem,”Science 292, 472-475 (2001).[9] T. Albash and D. A. Lidar, ”Adiabatic quantum computation,” Rev. Mod. Phys. 90, 015002(2018).[10] M. W. Johnson, et al., ”Quantum annealing with manufactured spins”, Nature 473, 194–198(2011).[11] R. D. Somma, S. Boixo, H. Barnum, and E. Knill, ”Quantum simulations of classicalannealing processes,” Phys. Rev. Lett. 101, 130504 (2008).[12] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, ”Quantum approximate optimizationalgorithm for MaxCut: A fermionic view,” Phys. Rev. A 97, 022304 (2018).[13] F. Gaitan and L. Clark, ”Ramsey numbers and adiabatic quantum computing,” Phys. Rev.Lett. 108, 010501 (2012).[14] F. Gaitan and L. Clark, ”Graph isomorphism and adiabatic quantum computing,” Phys.Rev. A 89, 022342 (2014).[15] P. L. McMahon, et al., ”A fully programmable 100-spin coherent Ising machine with all-to-all connections,” Science 354, 614-617 (2016).[16] K. Kudo, ”Constrained quantum annealing of graph coloring,” Phys. Rev. A 98, 022301(2018). 1217] A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: A Universal Concept inNonlinear Sciences (Cambridge University Press, 2003), Vol. 12.[18] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, ”Synchronization incomplex networks,” Phys. Rep. 469, 93-153 (2008).[19] M. A. Porter and J. P. Gleeson, “Dynamical systems on networks: A tutorial,”
Frontiersin Applied Dynamical Systems: Reviews and Tutorials (Springer, 2016), Vol. 4.[20] Z. Levnaji´c and A. Pikovsky, ”Phase Resetting of Collective Rhythm in Ensembles of Os-cillators,” Phys. Rev. E 82, 056202 (2010).[21] L. da F. Costa, et al., ”Analyzing and modeling real-world phenomena with complex net-works: a survey of applications,” Advances in Physics 60, 329-412 (2011).[22] V. Ja´cimovi´c and A. Crnki´c, ”Modelling mean fields in networks of coupled oscillators,” J.Geom. Phys. 124, 241-248 (2018).[23] A. Zakharova, “Chimera Patterns in Networks: Interplay between Dynamics, Structure,Noise, and Delay,” in
Understanding Complex Systems , edited by S. Kelso (Springer, 2020).[24] J.A. Acebr´on, L.L. Bonilla, C.J.P. Vicente, F. Ritort, and R. Spigler, ”The Kuramotomodel: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137(2005).[25] F.A. Rodrigues, T.K.D. Peron, P. Ji, and J. Kurths, ”The Kuramoto model in complexnetworks,” Phys. Rep. 610, 1-98 (2016).[26] A. Crnki´c and V. Ja´cimovi´c, ”Exploring complex networks by detecting collective dynamicsof Kuramoto oscillators,” in
Proceedings of the OPTIMA-2017 Conference , edited by Yu.G. Evtushenko, et al. (CEUR-WS, 2017), pp. 146-151.[27] C. Zankoc, D. Fanelli, F. Ginelli, and R. Livi, ”Desynchronization and pattern formationin a noisy feed-forward oscillator network,” Phys. Rev. E 99, 012303 (2019).[28] B. Pietras, N. Deschle, and A. Daffertshofer, ”First-order phase transitions in the Kuramotomodel with compact bimodal frequency distributions,” Phys. Rev. E 98, 062219 (2019).[29] Z. Levnaji´c, ”Emergent multistability and frustration in phase-repulsive networks of oscil-lators,” Phys. Rev. E 84, 016231 (2011).[30] Z. Levnaji´c, ”Evolutionary design of non-frustrated networks of phase-repulsive oscillators,”Scientific Reports 2, 967 (2012).[31] D. Goldstein, M. Giver, and B. Chakraborty, ”Synchronization patterns in geometricallyfrustrated rings of relaxation oscillators,” Chaos 25, 123109 (2015).[32] S. Astakhov, A. Gulai, N. Fujiwara, and J. Kurths, ”The role of asymmetrical and repulsivecoupling in the dynamics of two coupled van der Pol oscillators,” Chaos 26, 023102 (2016).[33] C. W. Wu, ”Graph Coloring via Synchronization of Coupled Oscillators”, IEEE Transac-tions on Circuits and Systems 45, 974-978, 1998.[34] J. Wu, L. Jiao, R. Li, and W. Chen, ”Clustering dynamics of nonlinear oscillator network:Application to graph coloring problem”, Physica D 240, 1972-1978, 2011.1335] A. V. Novikov and E. N. Benderskaya, ”Oscillatory Neural Networks Based on the Ku-ramoto Model for Cluster Analysis”, Pattern Recognition and Image Analysis 24, 365–371,2014.[36] A. Parihar, N. Shukla, M. Jerry, S. Datta, and A. Raychowdhury, ”Vertex coloring of graphsvia phase dynamics of coupled oscillatory networks”, Scientific Reports 7, 911, 2017.[37] G. Chartrand and P. Zhang,
Chromatic Graph Theory (Chapman and Hall/CRC, 2008).[38] T. R. Jensen, and B. Toft,
Graph coloring problems , Vol. 39, John Wiley & Sons (2011).[39] R. Ganguli, and R. Siddhartha, ”A study on course timetable scheduling using graph col-oring approach”, International Journal of Computational and Applied Mathematics 12,469-485 (2017).[40] K. I. Aardal, et al., ”Models and solution techniques for frequency assignment problems”,Annals of Operations Research 153, 79-129 (2007).[41] N. Gvozdenovi´c and M. Laurent, ”The operator Ψ for the chromatic number of a graph”,SIAM Journal on Optimization 19, 572-591 (2008).[42] J. Govorˇcin, N. Gvozdenovi´c, and J. Povh, ”New heuristics for the vertex coloring problembased on semidefinite programming”, Central European Journal of Operations Research 21,13-25 (2013).[43] A. Jabrayilov, and P. Mutzel, ”New integer linear programming models for the vertexcoloring problem”, In Latin American Symposium on Theoretical Informatics (pp. 640-652). Springer (2018).[44] D. E. Knuth, ”The sandwich theorem,” The Electronic Journal of Combinatorics 1, 1 (1994).[45] D. Eppstein, ”Small maximal independent sets and faster exact graph coloring”, J. GraphAlgorithms Appl. 7, 131-140 (2003).[46] E. Malaguti and P. Toth, ”A survey on vertex coloring problems”, International transactionsin operational research 17, 1-34 (2010).[47] R. Lewis,
A Guide to Graph Colouring: Algorithms and Applications (Springer, 2016).[48] R. K. Sundaram,
A First Course in Optimization Theory (Cambridge University Press,1996).[49] K. V. Mardia and P. E. Jupp,
Directional Statistics (John Wiley & Sons, 1999).[50] R. ˇSkrekovski, ”List improper colourings of planar graphs,” Combin. Probab. Comput. 8,293-299 (1999).[51] B. Mohar, ”Circular colorings of edge-weighted graphs,” J. Graph Theory 43, 107-116(2003).[52] A. Mishra, S. Banerjee, and W. Arbaugh, ”Weighted coloring based channel assignment forWLANs,” Mobile Comput. Commun. Rev. 9, 19-31 (2005).[53] U. Ascher and L. Petzold,
Computer Methods for Ordinary Differential Equations andDifferential-Algebraic Equations (SIAM, 1989).14
Full proofs for the weighted and non-weighted cases
In the appendices we report more complete and rigorous proofs for the statements made in themain text. We being by the simpler non-weighted case by noting that K -frustration along thelink ( i, j ) is defined as f Kij ( ϕ i , ϕ j ) = p K ( ϕ i − ϕ j ) , where p K is a trigonometric polynomial: p K ( x ) = 1 + c K ( K cos x + ( K −
1) cos 2 x + · · · + 2 cos( K − x + cos Kx ) , (A1)where the constants c K are chosen in such a way to have p K ( πK +1 ) = 0 satisfied.Notice several properties of polynomials p K :P1. For each K , p K is a trigonometric polynomial of the degree K .P2. p K are even functions.P3. p K are 2 π -periodic functions.P4. Functions p K are nonnegative, i.e. p K ( x ) ≥ x .P5. For each K the function p K has exactly K local minima on [0 , π ]. These local minimaare located at points πm +1 for m = 1 , . . . , K . (It is important to notice that there are no localminima at zero.) Definition.
We say that phases ϕ , . . . , ϕ N are in k -regular configuration for some k ≤ N ,if each ϕ i coincides with one of the points z + sπk , for some z ∈ [0 , π ] and s = 0 , . . . , k − Proposition 1.
1. For any graph Γ and for all integers K the functions F K are nonnegative.2. There exists an integer m ≤ N −
1, such that the value of F m at its global minimumequals zero.3. Suppose that F m ( ϕ , . . . , ϕ N ) = 0. Then, the configuration of phases ϕ , . . . , ϕ N is m + 1-regular. Proof.
1. Obviously, for each K the function F K is nonnegative, as it is the sum of nonnegativefunctions (A1), see the property P4.2. Due to property P5, polynomials p N − have precisely N local minima at points πm , m =1 , . . . , N . Let ϕ , ϕ = ϕ + πN , ϕ = ϕ + πN , . . . , ϕ N = ϕ N − + πN . Then, for each pair( i, j ), we have f N − i,j ) ( ϕ i , ϕ j ) = p N − ( ϕ i − ϕ j ) = 0. We have F ( ϕ , ϕ , . . . , ϕ N ) = 0. Hence, theglobal minimum of F N − is zero. In other words, the function F N − equals zero at N -regularconfiguration.3. Suppose that F m ( ϕ , . . . , ϕ N ) = 0. Then all terms in the function F m are equal to zero,i.e. for ∀ ( i, j ) ∈ E, f m ( i,j ) ( ϕ i , ϕ j ) = p m ( ϕ i − ϕ j ) = 0. Then, from the property P5 it follows that ϕ i − ϕ j = sπm +1 , for some s = 0 , , . . . , m . Proposition 2.
A proper coloring of a graph Γ with m colors exists if and only if the global minimum of F m − equals zero. Proof.
Suppose that for some ϕ , . . . , ϕ N , F m − ( ϕ , . . . , ϕ N ) = 0. Then, for each ( i, j ) ∈ E , we havethat f m − ( ϕ i , ϕ j ) = p m − ( ϕ i − ϕ j ) = 0. Now, from the property P5 we have that ϕ i − ϕ j = sπm for some s = 0 , , . . . , m −
1. This means that points e iϕ , . . . , e iϕ N are located at m pointson the unit circle, where an angle between any two consecutive points equals πm . Moreover, if15 i, j ) ∈ E , then ϕ i = ϕ j . Now, we assign a different color to each of m points on the unit circle.Then, if ( i, j ) ∈ E , the corresponding nodes q i and q j are colored in different colors.On the other side, suppose that graph Γ can be properly colored using m colors. Pick m points on the unit circle, so that an angle between any two consecutive is πm − and assign one of k colors to each point. To each node q i assign a value ϕ i , such that the point e iϕ i is colored into thesame color as q i . Suppose that ( i, j ) ∈ E . Then nodes q i and q j are colored into different colorsand, hence, ϕ i − ϕ j = sπm for some s = 1 , . . . , m −
1. Then, f m − i,j ) ( ϕ i , varphi j ) = 0. In otherwords, m − i, j ) ∈ E is zero, and, hence, the total m − F m − ( ϕ , . . . , ϕ N ) = 0.From Proposition 2 we immediately obtain the following Corollary.
The chromatic number of graph Γ is a minimal integer m , such that the global minimum of F m − equals zero.Finally, we extend above proofs to the weighted case. To do that, we prove the Propositionthat guarantees that our methods for minimization of functions F m − W will end up at m -regularconfigurations. This means that minimization method for F m − W will yield a certain coloring ofthe graph using m colors. Proposition 3.