Vertex coloring of graphs via phase dynamics of coupled oscillatory networks

Download PDF Article Open access Published: 19 April 2017 Vertex coloring of graphs via phase dynamics of coupled oscillatory networks Abhinav Parihar  ORCID: orcid.org/0000-0001-8203-58881,

Nikhil Shukla2, Matthew Jerry2, Suman Datta2 & …Arijit Raychowdhury1 Show authors Scientific Reports volume 7, Article number: 911 (2017) Cite this article

9081 Accesses

99 Citations

86 Altmetric

Metrics details

Subjects Electronic devicesMathematics and computing An Author Correction to this article was published on 12 April 2018

This article has been updated

While Boolean logic has been the backbone of digital information processing, there exist classes of computationally hard problems wherein this paradigm is fundamentally inefficient. Vertex

coloring of graphs, belonging to the class of combinatorial optimization, represents one such problem. It is well studied for its applications in data sciences, life sciences, social

sciences and technology, and hence, motivates alternate, more efficient non-Boolean pathways towards its solution. Here we demonstrate a coupled relaxation oscillator based dynamical system

that exploits insulator-metal transition in Vanadium Dioxide (VO2) to efficiently solve vertex coloring of graphs. Pairwise coupled VO2 oscillator circuits have been analyzed before for

basic computing operations, but using complex networks of VO2 oscillators, or any other oscillators, for more complex tasks have been challenging in theory as well as in experiments. The

proposed VO2 oscillator network harnesses the natural analogue between optimization problems and energy minimization processes in highly parallel, interconnected dynamical systems to

approximate optimal coloring of graphs. We further indicate a fundamental connection between spectral properties of linear dynamical systems and spectral algorithms for graph coloring. Our

work not only elucidates a physics-based computing approach but also presents tantalizing opportunities for building customized analog co-processors for solving hard problems

efficiently.

photonics-inspired clock models Article Open access 29 April 2022 An Ising Hamiltonian solver based on coupled stochastic phase-transition nano-oscillators Article 26 July 2021 Introduction

The semiconductor industry is pivoted upon the Von Neumann computer architecture which implements the “Turing Machine” model of computation with a clear distinction between processing units

and memory. Computation is carried out through a sequence of instructions with periodic loads and stores to the memory. On the contrary, computation in nature, our brain included, follows a

radically different approach. Processing is distributed in all parts of the machine; memory and processors are integrated; clear distinguishable atomic instructions are replaced by

continuous time dynamics; and information is encoded in physically meaningful quantities instead of their symbolic interpretations. In spite of the success of the von Neumann computing

architecture, its limitations become apparent1,2 when dealing with certain classes of problems such as associative computing3,4, optimizations5, pattern matching and recognition6. This has

motivated active research in alternative computing models, where dynamical systems have been shown to provide a fundamentally new platform to address these increasingly important problem

classes7,8,9,10,11.

In this paper, we report experimental evidence and the corresponding theoretical foundation for harnessing the continuous time dynamics of a system of coupled relaxation oscillators to solve

vertex coloring of a random graphs, a combinatorial optimization problem of large-scale importance. Combinatorial optimizations represent a problem class where an optimal value of a

function, or its optimal point, needs to computed within a domain set which is discrete or combinatorial. Vertex coloring of graphs is a combinatorial optimization problem which is NP-hard

(non-deterministic polynomial-time hard), unless P = NP. This means that the best algorithms end up searching the whole domain set for at least some problem instances. Vertex coloring is

also one of the most studied NP-hard combinatorial optimization problems not only for its significance in computational theory but also for its many real world applications like fault

diagnosis12, scheduling13,14,15, resource allocation16; Such problems are believed to be solved, or approximated, efficiently in natural processes because they can explore the solution space

in a massively parallel manner17. In fact, for any deterministic system to be able to solve such hard problems, be it a sequential deterministic Turing machine or a continuous time

dynamical system, exponential resources are required which can be in terms of time, hardware components, maximum magnitudes of variables or their precision7,18,19,20.

In the last three decades, dynamical systems as well as hardware implementations have been proposed to solve NP problems, many of which implement some form of algorithm for solving, or

approximating, such problems. Important attempts include Quantum computers21, Cellular automata22, Hopfield networks5, Ising model formulations23, chaotic nonlinear attractor systems7,

iterated projections24, Memcomputing25 and stochastic searching using non-repeating phase relations among oscillators26. All these approaches, except iterated maps, are based on the idea of

interconnected “nodes” which exchange, store and process information among themselves. One particular variant of such dynamical systems is based on a network of coupled oscillators whose

phase and frequency dynamics can be exploited to encode system states3,27,28. These computational kernels, which have been claimed to mimic the spiking networks of the human brain, have been

successfully used in associative computing, demonstrating significant improvements in energy-efficiency and performance when applied to video analytics29. Theoretical models of coupled

sinusoidal oscillators, (Kuramoto models30,31), or van der pol oscillators32,33 have often been used to theoretically study asymptotic limits; but their physical implementations and

experimental evidence of such networks to solve computationally hard problems have remained elusive.

Here, we establish that a system of coupled relaxation oscillators fabricated using Vanadium dioxide (VO2) metal-insulator-transition devices and coupled capacitively, can lead to system

dynamics on which vertex coloring of unweighted and undirected graphs (hitherto referred to as the graph coloring problem) can be successfully mapped (Fig. 1a). A VO2 oscillator consists of

a VO2 device connected in series with a conductance g s (with a loading capacitor c L in parallel) and the output node is the node between the VO2 device and the conductance (Fig. 1a). Such

a simple series circuit shows self-sustained relaxation oscillations. In previous works by authors, properties of two such coupled Vanadium Dioxide (VO2) based relaxation oscillators have

been analyzed28,34 and pairwise coupled circuits have been proposed for basic computation tasks3,35, but using complex network connections with VO2 oscillators (or any other kind of

oscillators) for more complicated computing tasks have been challenging in theory as well as in experiments. We demonstrate experimentally and using simulations that when such relaxation

oscillators are coupled using only capacitances in a manner topologically equivalent to an input graph, their steady state phases can be used to approximate the solution of the NP-hard

minimum graph coloring problem. For this, we propose a reformulation of the graph coloring problem where instead of finding a color assignment for each node, the objective is to find a

circular ordering or circular permutation of the nodes such that the same colored nodes appear together in the ordering. Such a reformulation preserves the hardness of the problem and is

useful for interpreting the output of our circuit (Fig. 1a). We show analytically that the dynamics of such a coupled relaxation oscillator system is intrinsically connected to spectral

algorithms for graph coloring36,37,38, which use eigenvectors of adjacency matrix of the input graph to approximate the solutions. Alternatively, the permutation of steady state phases of

coupled relaxation oscillators depends on eigenvectors of the adjacency matrix in the same way as have been used by spectral algorithms for graph coloring (Fig. 1b). A programmable circuit

for a such a coupled oscillator system, where the oscillators are coupled in a graph with adjacency matrix, A and coupling capacitance c c is shown in Fig. 2. Such an all-to-all connected

circuit can be laid out in 2-dimensions with each connection controlled by a switch (transistor). The gate voltages of transistors can be used to provide the input, viz. the adjacency matrix

A of the problem graph instance. Our simulation results show that the hardness of problem instances has, on average, expected effects on important metrics of solutions found using such a

circuit like the number of colors detected and the settling time.

Overview of the circuit and system dynamics. (a) Overview of the proposed system for vertex coloring and a simulation example. First step is a coupled relaxation oscillator circuit where the

oscillators are composed of a series combination of VO2 device and a resistor (with a loading capacitor in parallel), and are connected in a graph using capacitors. The equivalent circuit

diagram of the VO2 oscillator is shown using an internal capacitance c i and a phase changing conductance g(m/i) which switches between metallic conductance g m and insulating conductance g

i . An example 3-partite graph is simulated and the relative phases of these oscillators are shown in a phase diagram which shows vertex color-sorting in phase, and can be used to calculate

vertex-coloring with O(n2) complexity. (b) The circuit is composed of VO2 oscillators capacitively coupled in a network same as the input graph. The final order of phases, or charging

spikes, of the oscillators is related to the eigenvectors of the adjacency matrix of the input graph which in turn are related to the solution of the graph coloring problem.

imageFigure 2

Coupled oscillator circuit schematic. A circuit of 4 coupled oscillators with capacitive connections between oscillators controlled using switches corresponding to the adjacency matrix A and

coupling capacitance c c . The subscripts denote the corresponding entries in A. Note that A ij  = A ji , A ii  = 0 and Aij ∈ {0,1}.

It is well known that eigen properties of the coefficient matrix in the evaluation equation of a dynamical system determine important structural properties of the system including stability,

bifurcation, energy minima(s) and overall system dynamics. In this report, we provide a theoretical bridge and experimental evidence that a dynamical system whose coefficient matrix

inherits properties of the incidence matrix of a graph, can indeed emulate spectral graph algorithms just through its time evolution. We envision such dynamical systems to provide

foundational paradigms in the development of next-generation computational accelerators and kernels.

The objective of graph coloring or vertex coloring is to assign one color (out of total k colors) to each vertex of an undirected graph such that no two adjacent vertices receive the same

color. A graph coloring that minimizes the number of colors k is called minimum graph coloring. The minimum k for which a correct coloring is possible is called the chromatic number of the

graph. A graph which can be colored using at most k colors is called a k-partite graph. A k-partition of a set, like the set of nodes, is a grouping of the elements of the set into k groups.

Hence, a vertex coloring with k colors is a k-partition. We reformulate the objective of finding a color assignment to finding a circular permutation of nodes such that the same colored

nodes appear together. We refer to this reformulation as vertex color-sorting and the corresponding optimal version as minimum vertex color-sorting. Calculating a color assignment from a

color-sorting is O(n2). This is because if P is the permutation matrix for the color-sorting and A is the symmetric adjacency matrix of the graph, then the color assignments can be found by

observing the ‘0’ diagonal blocks in the matrix PAPT. This makes vertex color-sorting as hard as vertex coloring (see Supplementary Section 4). Using this method, any permutation P gives a

correct color assignment but a better permutation gives lesser number of colors, and an optimal color-sorting permutation gives the minimum number of colors. In the proposed coupled

oscillator system, each oscillator represents a vertex (or node) of the graph. Any two nodes connected in the original graph by an edge (as indicated by a ‘1’ in A), are capacitively coupled

in the hardware implementation. As the coupled system evolves, the relative phases of the oscillators are ordered, and we observe that the relative ordering of the phases approximates

minimum vertex-color-sorting of the original graph.

We construct a relaxation oscillator by exploiting the electrically induced large and abrupt change in resistance across the insulator-to-metal transition (IMT) in Vanadium Dioxide

(VO2)39,40, and stabilizing it with a negative feedback from a series conductance g s (details of the single oscillator dynamics have been elucidated in our previous work28,34). Since the

IMT in VO2 is a materials-level manifestation of hysteretic, resistive threshold switching behavior41,42 critical to realizing relaxation oscillatory action43, VO2 based oscillators present

a compact, scalable, and potentially low power solution35 to realizing the fundamental building block of the graph coloring hardware. Further, we use a non-dissipative capacitive coupling

scheme to connect the oscillators and achieve frequency synchronization.

Figure 3 shows two representative configurations of graphs (Fig. 3a: delta configuration; Fig. 3b: cross-connected ring configuration) along with their equivalent implementations using

coupled oscillators. The respective time domain waveforms of the oscillators (Fig. 3c,d) reveal a unique relationship among the phases of the oscillators: there is a distinct non-zero phase

difference between any two directly coupled oscillators. This is because the nature of capacitive coupling among the relaxation oscillators ensures that two adjacently connected oscillators

will tend to force each other out of phase; a rigorous mathematical treatment of the phase dynamics between two coupled VO2 based oscillators has been detailed in an earlier work34.

Additionally, when an oscillator is connected to multiple other nodes, the net phase of the oscillator is the aggregate of the ‘repelling effect’ of all the other connected oscillators. As

such, in the light of vertex coloring, such a circuit is expected to result in output phases which are clustered by color, i.e. oscillators with the same color have phases which are close

together. But such an interpretation of oscillator phases is weak and is difficult to apply in most cases of graphs where the corresponding circuit output either does not have well clustered

phases or have incorrect clusters. Our interpretation of outputs as color-sorting solves all these problems and is well defined. As will be discussed in the next section, the combined

repelling effect in a network of oscillators gives special properties to the order of phases of oscillators in steady state, viz. they approximate minimum vertex-color-sorting. As discussed

earlier, this steady state ordering of phases is then used to calculate a vertex-coloring.

Phase dynamics of synchronized VO2 based capacitively coupled relaxation oscillators. (a,b) Schematics of two representative configurations (a: delta configuration; b: cross-connected ring)

of capacitively coupled VO2 based oscillators, and their corresponding graphs. (c,d) Time domain waveforms from experiment and simulations for the two coupled oscillator configurations in

(a,b), respectively, showing that while the oscillators are synchronized in frequency, no two directly coupled oscillators are in-phase. This important property of the coupled oscillator

system enables graph coloring. (e,f) Time averaged XOR of thresholded outputs of oscillators (each w.r.t. oscillator number 1), and respective polar phase plots showing steady state relative

phases detected using PFDs. The XOR values are normalized with respect to the maximum value.

Since the oscillators in this work are non-sinusoidal in nature, the steady state phase differences among the coupled oscillators can be calculated using phase-frequency detectors (PFDs) or

the time-averaged XOR metric35. The time-averaged XOR measure of any two oscillator outputs is calculated by first thresholding the outputs to binary valued waveforms and then taking the

time average of absolute difference of these thresholded waveforms over the complete steady-state periodic orbit of the system. The time-averaged XOR metric is proportional to the absolute

value of phase difference between the oscillators, and hence it does not differentiate between lead or lag. Figure 3e,f show the relative phases detected using a PFD (shown using the polar

phase plots) and the XOR measures of each oscillator with respect to a common reference oscillator (shown as bar graphs).

Next, we experimentally investigate the coloring of some other graph configurations with up to five vertices, using the system of VO2 based oscillators (Fig. 4). The coupled oscillators are

configured to represent the respective graphs as discussed earlier, and the corresponding values of the time-averaged XOR along with the respective phase plots of the oscillators are shown

in Fig. 4. It can be observed that the hardware is able to optimally color all the graphs investigated here.

Experimental results of graph coloring using coupled VO2 based oscillators. Various graph configurations, and their experimentally obtained solutions (PFD outputs and XOR values) using the

coupled relaxation oscillator system. After mapping the graphs onto the coupled oscillator hardware, the steady state order of phases of oscillators is used as a color-sorting and a color

assignment is calculated.

The mathematical model of the circuit is created taking into consideration the fact that the fundamental switching time constant for the phase transition (i.e., IMT and MIT) in correlated

materials like Vanadium Dioxide (VO2) occurs on extremely fast time-scales (~75 fs)44,45. Even for electrically induced phase transitions, the expected time constants in scaled VO2 devices46

are a few orders of magnitude smaller than the RC time constants observed in the oscillator circuit. As such, the phase transition is modeled as an abrupt instantaneous switching. The VO2

devices switch between a low resistance metallic state with conductance g m and a high resistance insulating state with conductance g i based on the voltage v d across their two terminals.

On increasing v d the device switches to a metallic state (insulator-to-metal (IMT) transition) after a threshold v h , and on decreasing v d below v l the device switches back to an

insulating state (metal-to-insulator (MIT) transition). Here v h  > v l and v h  − v l defines the hysteresis in switching. Consider a supply voltage which is applied across the series

combination of such a hysteretic device and a conductance g s where the subscript s denote a series conductance. Without loss of generality, we assume that at t = 0 the device is in high

resistance state and the voltage drop across the device v d  = 0. The internal capacitance of the device charges up and v d increases and eventually crosses the threshold v h . Due to this

the device transitions into a metallic state which causes the internal capacitance of the device to discharge and reduces v d which finally drops below v l . This causes the device to switch

back to the insulating state resulting in oscillations with piecewise linear dynamics. In case of the coupled oscillator circuit, a loading capacitance c L of appropriate magnitude is

required as shown in Fig. 1a for correct circuit operation.

The dynamics of a single oscillator (Fig. 1a) can be written as the following piecewise differential equation:

where c is the lumped capacitance of device along with the loading capacitance and parasitics, s ∈ {0,1} is the state of system - charging (denoted by 1) or discharging (denoted by 0), and

g(s) is the net path conductance in state s, with g(s) = g s  + g i s. If the voltage v is normalized to v dd then p(s) = s. The dynamics of a circuit of identical coupled relaxation

oscillators can be described using the following matrix differential equation:

}}(t)={({C}_{i}+{C}_{c}+{C}_{l})}^{-1}[-G({\boldsymbol{s}}){\boldsymbol{x}}({\boldsymbol{t}})+{g}_{i}{\boldsymbol{s}}]$$ (2)

Here, x is the vector of all voltages (normalized to v dd ), C i is a diagonal matrix with the diagonal elements equal to the internal capacitances of the corresponding oscillator nodes, C c

is the coupling capacitance matrix with diagonal elements equal to the sum of all the coupling capacitances connected to the corresponding nodes and off-diagonal elements equal to the

coupling capacitances of corresponding pair of nodes with negative sign, s ∈ {0,1}n is the vector of states of all oscillators, G(s) is a state dependent diagonal matrix with

\(diag(G({\boldsymbol{s}}))={g}_{s}+{g}_{i}{\boldsymbol{s}}\), and C l is a diagonal matrix corresponding to the extra loading capacitors. These loading capacitors effectively add to the

internal capacitance and are chosen such that diag(C c  + C l ) is constant. When all oscillators have equal internal capacitances c i and equal coupling capacitances c c , then C i  = c i I

where I is the identity matrix, C c  = c c L, L being the laplacian matrix of the graph with L = D − A where D is a diagonal matrix of degrees of vertices and A is the adjacency matrix of

the graph. In such a case, a simple choice of C l is C l  = c c (nI − D). These oscillators have very high charging rate (due to the high conductance of the metallic state of the VO2

devices) and low discharging rates (due to relatively high resistance of the pull-down resistors), and as such, the phases of oscillators and their permutation can be read by observing the

relative positions of the charging spikes.

Considering the system of (2), if there exists a limit cycle where the system settles to a certain order of charging spikes, then the order of charging spikes is same as the order of

components of the state vector x in state s = 0 (Fig. 5a), where the orders are considered unique up to a circular permutation and 0 is a vector of zeros. In other words, orders which are

circular permutations of each other are considered same. The system of (2) in state s = 0 is a linear dynamical system with all negative eigenvalues and the asymptotic order of components of

the state vector x is determined by the asymptotic direction of the system trajectory. A representative figure of system trajectories in such a linear dynamical system in two dimensions is

shown in Fig. 5b in which e1 and e2 are the eigenvectors with distinct negative eigenvalues of the coefficient matrix B of a linear dynamical system

\(\dot{{\boldsymbol{x}}}=B{\boldsymbol{x}}\), where x = {x,y}. The eigenspaces E1 and E2 are the lines along e1 and e2 respectively and PE1 and PE2 are projection matrices for these

eigenspaces. Considering (2), the eigenvectors of coefficient matrix (C i  + C c  + C l )−1 with the least negative eigenvalues are, in fact, same as the eigenvectors of adjacency matrix A

with most negative eigenvalues (Supplementary Proposition 1). These eigenvectors determine the asymptotic order of components of the state vector x in the following way. Let T(x) represent

the order of components of vector x. Referring to Fig. 5b, let x0 be the initial starting point of the system, E1 the eigenspace with least negative eigenvalue, E2 with the next higher

eigenvalue and so on. Then the asymptotic order of components is same as the order of components of PE1x0. In case PE1x0 has some components which are equal (with respect to Fig. 5b, it

means if e1 lies along the x = y line) then the order among those equal components is decided by the order of PE2x0 and so on (see Supplementary Section 2 for analytical derivations). If

this operation of combining two orders is represented by \(\oplus ^{\prime} \) where the first order is preferred over the second, then the asymptotic order can be written

as:

System dynamics and asymptotic permutation. (a) Simulation waveforms of a circuit connected in the 3-partite graph shown in Fig. 1a. The gray regions show the time when the system is in

state s = 0. (b) Representative figure showing the relation between the asymptotic order of components of the state vector and the eigenspaces in a two-dimensional linear system where the

coefficient matrix has negative eigenvalues.

The eigenvectors with least negative eigenvalues determine not just the asymptotic order in state s = 0 but also the limit cycle. To obtain an intuitive understanding of how the system

settles to a limit cycle with the correct color-sorting, i.e. steady state oscillations where order of charging spikes become constant and equal to a correct color-sorting, we note that such

a limit cycle exists if:

The state s = 0 does not change T(x) when T(x) is a correct color-sorting.

The charging states s ≠ 0 change T(x) by only a circular permutation.

The state s = 0 does not change T(x) even when T(x) is a circular permutation of the correct color-sorting for which property 1 holds.

Property 1 is true when the system is at state x such that T(x) corresponds to the order determined by the least negative eigenvalues (from (3)) because then it is same as the asymptotic

order of components and hence does not change. Also, the eigenvectors with the least negative eigenvalues have a property, as used in spectral algorithms for graph coloring36,37,38, that

their components corresponding to same color tend to be equal (or close) for dense graphs, and these components disperse as graphs become sparse. For instance, in a complete 3-partite graph

with arbitrary partition, the eigenvectors have exactly equal components for nodes belonging to the same partition subset (color class). Hence, the order determined by these least negative

eigenvectors (or eigenspaces) as given by (3) will correspond to a correct color-sorting with minimum number of colors for dense graphs and the number of colors would increase for sparser

graphs. This can also be seen in the light of perturbation theory of matrices. Any sparse k-partite graph (or equivalently its adjacency matrix) can be obtained from a complete k-partite

graph with the same partition structure by removing some edges (or changing ‘1’s in its adjacency matrix to ‘0’s). Such a “perturbation” of adjacency matrix has an equivalent effect of

rotating its eigenvectors by some angles47 which depend on the magnitude of perturbation, which in this case is the number of edges removed. We empirically find properties 1–3 to be true for

complete k-partite graphs which are the densest k-partite graphs (see Supplementary Section 5 for an analytical discussion), and considering sparse graphs as perturbations of complete

partite graphs, we can say that vertex color-sorting using the coupled relaxation oscillator circuit becomes less optimal as graphs become sparse. This is known to be true for spectral and

other coloring algorithms as well that k-colorable dense graphs are easier to color than sparser ones.

We simulate the dynamical system as described by (2) for random graph instances of 3-colorable graphs. The initial conditions are chosen at random \({{\boldsymbol{x}}}_{0}\in

{[0.35,0.65]}^{n}\) and s = 1 at t = 0 because all oscillators are in charging state when the power is switched on. Without loss of generality, v l and v h are chosen as 0.2 V and 0.8 V

respectively with a supply voltage of 1 V. We use a random graph generation model G(n,k,i) to generate instances of colorable k-partite graphs with total n nodes. The graphs are generated by

first choosing a random k-partition of n nodes, then creating a complete k-partite graph with this k-partition and finally removing random i number of edges from this complete graph.

Average connectivity is defined as the ratio of total number of edges in the generated graph G(n,k,i) to the total number of edges in the complete k-partite graph with the same partition.

As is true with hard problems, even in graph coloring problems no heuristic algorithm works best for all graph instances48. Also different heuristics work better for different instances, and

hence no single order parameter can account for the hardness of an instance of a graph coloring problem49,50,51. The most commonly used order parameter is average connectivity52. We use

this parameter to account for the hardness of the problems being solved and observe how a coupled relaxation oscillator network behaves for problems with varying levels of average

connectivity. Observations are made particularly about the cluster diameter, the number of colors detected and the settling time, which are defined as follows.

When the coupled oscillator circuit settles to a correct color-sorting, the phases of oscillators or nodes with the same color form a cluster for many graphs, esp. the dense graphs. The

maximum phase difference of two oscillators in the same cluster, i.e. with the same color, is called the cluster diameter. The number of colors detected is calculated using the order of

charging spikes at the end of the finite time period for which the circuit is simulated. Settling time is the defined as the time after which the number of colors detected does not change

till the end of the simulation time.

Figure 6 gives a visualization to how the order of charging spikes evolves with time for 3 different graphs of 20 nodes with decreasing average connectivity. All three graphs are 3-partite

with partition (8,2,10). For a single simulation instance, we note the order of charging spikes at various time instances and associate a unique number (within a simulation instance) to each

permutation. A plot of this permutation number with time shows how the order, or permutation, of the charging spikes evolves with time. Figure 6 shows this plot along with plot of the

number of colors detected using the order of charging spikes at various times. Figure 6a shows the typical case of a complete partite graph where the order of charging spikes settles quickly

to a correct color-sorting, and the number of colors detected falls quickly to the minimum number of colors (3 in this case). Figure 5b and c show graphs with lower connectivity but the

same partition structure. We make two observations. Firstly, even after the number of colors detected settles down, the permutation or order of charging spikes can evolve. Secondly, Fig. 6c

shows lower number of colors detected than Fig. 6b but the settling time is higher for Fig. 6c. As such, both settling time and number of colors detected can be considered as imperfect order

parameters for hardness of graph coloring just like average connectivity.

Effect of hardness on circuit behavior. (Above) Graph diagram and (below) waveforms of permutation number (a unique number corresponding to each permutation) of charging spikes and number of

colors detected with time. All graphs are 3-partite with partition (8,2,10) with varying levels of average connectivity: (a) 1.0, (b) 0.57 and (c) 0.48. All triangles in the graph are shown

with red edges and the rest edges are shown in blue.

Figure 7 shows the performance of such network on random graph instances. We generate 3-partite graphs using G(10,3,i), G(20,3,i) and G(30,3,i) with increasing values of i. The graphs which

become bipartite after removing i edges are discarded. Each graph instance is used 5 times with a new random initial condition for the relaxation oscillator circuit. Various metrics to

evaluate the circuit output are plotted against average connectivity. We see that as graphs become sparse with decreasing average connectivity, the cluster diameter increases (Fig. 7a). This

comparison is made only among those graphs where the final phases are clustered correctly into 3 clusters, i.e. the number of colors detected is 3. When graphs are dense and closer to being

complete partite, the possibility of them being optimally colored with 3 colors is high, and the settling time on average is less (Fig. 7b,c). As graphs become sparser, the number of colors

detected (Fig. 7b) as well as settling time (Fig. 7c) increase statistically on average. It also follows our intuitive understanding that hard computational problems remain hard even under

domain transformation, albeit with potential practical implications such as increased energy-efficiency and performance benefits of continuous time systems over their digital counterparts. A

comparison of colors detected from simulating the coupled oscillator network with that using Brelaz Heuristics53 (Fig. 7d) shows the effectiveness of the circuit as a tool to approximate

the minimum graph coloring problem. Number of colors detected by simulating sample graphs from the second DIMACS implementation challenge54 are shown in Table 1, where for certain instances

we note that the dynamical system outperforms heuristic algorithms.

Simulation results on random graph instances. (a) Maximum cluster diameter for those graphs for which 3 colors were detected. (b) (Top) Number of colors detected plotted against the average

connectivity. (Below) Mean colors detected in connectivity intervals. (c) (Top) Settling time plotted against the average connectivity. (Below) Mean settling time in connectivity intervals.

(d) Number of colors detected using the relaxation oscillator circuit plotted against number of colors detected using Brelaz heuristics for the random graph instances used in b and c. Each

graph instance represents a point whose coordinates are denoted by the pair of colors detected using the two methods. The size of dots represents the percentage of instances which lie at the

center of corresponding dot.

In this article, we have established that a system of capacitively coupled relaxation oscillators can perform graph coloring, which is a commonly studied and practically useful combinatorial

optimization problem. Further, the connection between system dynamics and the order of steady state phases of oscillators with spectral techniques for graph coloring has been discussed and

it shows an innate, yet, hitherto unexplored, connection between the time evolution of dynamical systems and computationally hard problems that have solutions or approximations in the

spectral domains.

The VO2 films have a thickness of 10 nm, and are epitaxially grown on (001) TiO2 using reactive oxide molecular beam epitaxy. The epitaxial mismatch between VO2 and TiO2 results in a tensile

biaxial strain of −0.9%.

The electrodes are patterned using contact lithography followed by electron beam evaporation of Pd/Au (20 nm/80 nm) and lift-off in RemoverPG at 70 °C. Next, the channel width and isolation

are defined by electron beam lithography followed by a CF4 dry etch. Finally, the resist is stripped with RemoverPG at 70 °C.

The oscillator circuits were simulated in Mathematica 10.2 for a finite time (1000 time units). Simulations were performed using default settings for NDSolve routines. The metal-insulator

transition events were detected using the inbuilt Mathematica event detection in NDSolve routines with default settings. For Brelaz heuristics, Mathematica routine for Brelaz heurstics was

used.

A correction to this article has been published and is linked from the HTML and PDF versions of this paper. The error has not been fixed in the paper.

References Theis, T. N. & Solomon, P. M. In Quest of the ‘Next Switch’: Prospects for Greatly Reduced Power Dissipation in a Successor to the Silicon Field-Effect Transistor. Proc. IEEE 98,

2005–2014, doi:10.1109/JPROC.2010.2066531 (2010).

Article  Google Scholar 

Zhirnov, V. V., Cavin, R. K., Hutchby, J. A. & Bourianoff, G. I. Limits to binary logic switch scaling - a gedanken model. Proc. IEEE 91, 1934–1939, doi:10.1109/JPROC.2003.818324 (2003).

Article  Google Scholar 

Shukla, N. et al. Pairwise coupled hybrid vanadium dioxide-MOSFET (HVFET) oscillators for non-boolean associative computing. In Electron Devices Meeting (IEDM), 2014 IEEE International

28.7.1–28.7.4, 10.1109/IEDM.2014.7047129 (2014).

Hoppensteadt, F. C. & Izhikevich, E. M. Synchronization of laser oscillators, associative memory, and optical neurocomputing. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.

62, 4010–4013, doi:10.1103/PhysRevE.62.4010 (2000).

CAS  Google Scholar 

Hopfield, J. J. & Tank, D. W. ‘Neural’ computation of decisions in optimization problems. Biol. Cybern. 52, 141–152 (1985).

CAS  PubMed  MATH  Google Scholar 

Chua, L. O. & Yang, L. Cellular neural networks: applications. IEEE Trans. Circuits Syst. 35, 1273–1290, doi:10.1109/31.7601 (1988).

Article  MathSciNet  Google Scholar 

Ercsey-Ravasz, M. & Toroczkai, Z. Optimization hardness as transient chaos in an analog approach to constraint satisfaction. Nat. Phys. 7, 966–970, doi:10.1038/nphys2105 (2011).

Article  CAS  Google Scholar 

Nikonov, D. E. et al. Coupled-Oscillator Associative Memory Array Operation for Pattern Recognition. IEEE J. Explor. Solid-State Comput. Devices Circuits 1, 85–93,

doi:10.1109/JXCDC.2015.2504049 (2015).

Article  Google Scholar 

Welser, J. J., Bourianoff, G. I., Zhirnov, V. V. & Cavin, R. K. The quest for the next information processing technology. J. Nanoparticle Res. 10, 1–10, doi:10.1007/s11051-007-9305-8 (2007).

Article  ADS  Google Scholar 

Srinivasan, G., Sengupta, A. & Roy, K. Magnetic Tunnel Junction Based Long-Term Short-Term Stochastic Synapse for a Spiking Neural Network with On-Chip STDP Learning. Sci. Rep. 6, 29545,

doi:10.1038/srep29545 (2016).

Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

Sharad, M., Augustine, C. & Roy, K. Boolean and non-Boolean computation with spin devices. In 11.6.1–11.6.4, doi:10.1109/IEDM.2012.6479026 (IEEE, 2012).

Akers, S. B. Fault Diagnosis as a Graph Coloring Problem. IEEE Trans. Comput. 23, 706–713, doi:10.1109/T-C.1974.224018 (1974).

Article  MATH  Google Scholar 

Leighton, F. T. A graph coloring algorithm for large scheduling problems. J. Res. Natl. Bur. Stand. 84, 489–506, doi:10.6028/jres.084.024 (1979).

Article  MathSciNet  MATH  Google Scholar 

Zufferey, N., Amstutz, P. & Giaccari, P. Graph colouring approaches for a satellite range scheduling problem. J. Sched. 11, 263–277, doi:10.1007/s10951-008-0066-8 (2008).

Article  MathSciNet  MATH  Google Scholar 

Xizheng, Z. & Yaonan, W. New mixed broadcast scheduling approach using neural networks and graph coloring in wireless sensor network. J. Syst. Eng. Electron. 20, 185–191 (2009).

Google Scholar 

Woo, T.-K., Su, S. Y. W. & Newman-Wolfe, R. Resource allocation in a dynamically partitionable bus network using a graph coloring algorithm. IEEE Trans. Commun. 39, 1794–1801,

doi:10.1109/26.120165 (1991).

Article  Google Scholar 

Yan, L. et al. Some massively parallel algorithms from nature. Wuhan Univ. J. Nat. Sci. 7, 37–46.

Vergis, A., Steiglitz, K. & Dickinson, B. The complexity of analog computation. Math. Comput. Simul 28, 91–113, doi:10.1016/0378-4754(86)90105-9 (1986).

Article  MATH  Google Scholar 

Siegelmann, H. T. Computation Beyond the Turing Limit. Science 268, 545–548, doi:10.1126/science.268.5210.545 (1995).

Article  ADS  CAS  PubMed  Google Scholar 

Ercsey-Ravasz, M. & Toroczkai, Z. The Chaos Within Sudoku. Sci. Rep. 2, doi:10.1038/srep00725 (2012).

Shor, P. Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer. SIAM J. Comput. 26, 1484–1509, doi:10.1137/S0097539795293172 (1997).

Article  MathSciNet  MATH  Google Scholar 

Theory and applications of cellular automata: including selected papers, 1983–1986 (World Scientific, 1986).

Lucas, A. Ising formulations of many NP problems. Interdiscip. Phys. 2, 5, doi:10.3389/fphy.2014.00005 (2014).

Google Scholar 

Elser, V., Rankenburg, I. & Thibault, P. Searching with iterated maps. Proc. Natl. Acad. Sci. 104, 418–423, doi:10.1073/pnas.0606359104 (2007).

Article  ADS  MathSciNet  CAS  PubMed  PubMed Central  MATH  Google Scholar 

Traversa, F. L., Ramella, C., Bonani, F. & Di Ventra, M. Memcomputing NP-complete problems in polynomial time using polynomial resources and collective states. Sci. Adv. 1,

e1500031–e1500031, doi:10.1126/sciadv.1500031 (2015).

Article  ADS  PubMed  PubMed Central  Google Scholar 

Mostafa, H., Müller, L. K. & Indiveri, G. An event-based architecture for solving constraint satisfaction problems. Nat. Commun. 6, 8941, doi:10.1038/ncomms9941 (2015).

Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

Izhikevich, E. M. Weakly pulse-coupled oscillators, FM interactions, synchronization, and oscillatory associative memory. IEEE Trans. Neural Netw. 10, 508–526, doi:10.1109/72.761708 (1999).

Article  CAS  PubMed  Google Scholar 

Shukla, N. et al. Synchronized charge oscillations in correlated electron systems. Sci. Rep. 4, doi:10.1038/srep04964 (2014).

Narayanan, V. et al. Video analytics using beyond CMOS devices. In 2014 Design, Automation Test in Europe Conference Exhibition (DATE) 1–5, doi:10.7873/DATE.2014.357 (2014).

Strogatz, S. H. From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Phys. Nonlinear Phenom. 143, 1–20, doi:10.1016/S0167-2789(00)00094-4

(2000).

Article  ADS  MathSciNet  MATH  Google Scholar 

Wu, J., Jiao, L., Li, R. & Chen, W. Clustering dynamics of nonlinear oscillator network: Application to graph coloring problem. Phys. Nonlinear Phenom. 240, 1972–1978,

doi:10.1016/j.physd.2011.09.010 (2011).

Article  ADS  MathSciNet  MATH  Google Scholar 

Kopell, N. & Somers, D. Anti-phase solutions in relaxation oscillators coupled through excitatory interactions. J. Math. Biol. 33, 261–280, doi:10.1007/BF00169564 (1995).

Article  MathSciNet  CAS  PubMed  MATH  Google Scholar 

Wu, C. W. Graph coloring via synchronization of coupled oscillators. IEEE Trans. Circuits Syst. Fundam. Theory Appl. 45, 974–978, doi:10.1109/81.721263 (1998).

Article  Google Scholar 

Parihar, A., Shukla, N., Datta, S. & Raychowdhury, A. Synchronization of pairwise-coupled, identical, relaxation oscillators based on metal-insulator phase transition devices: A Model Study.

ArXiv14082582 Nlin (2014).

Parihar, A., Shukla, N., Datta, S. & Raychowdhury, A. Exploiting Synchronization Properties of Correlated Electron Devices in a Non-Boolean Computing Fabric for Template Matching. IEEE J.

Emerg. Sel. Top. Circuits Syst. PP, 1–10, doi:10.1109/JETCAS.2014.2361069 (2014).

Aspvall, B. & Gilbert, J. Graph Coloring Using Eigenvalue Decomposition. SIAM J. Algebr. Discrete. Methods 5, 526–538 (1984).

MATH  Google Scholar 

Alon, N. & Kahale, N. A Spectral Technique for Coloring Random 3-Colorable Graphs. SIAM J. Comput. 26, 1733–1748, doi:10.1137/S0097539794270248 (1997).

Article  MathSciNet  MATH  Google Scholar 

McSherry, F. Spectral partitioning of random graphs. In 42nd IEEE Symposium on Foundations of Computer Science, 2001. Proceedings 529–537, 10.1109/SFCS.2001.959929 (2001).

Kim, H.-T. et al. Mechanism and observation of Mott transition in VO2 based two- and three-terminal devices. New J. Phys. 6, 52, doi:10.1088/1367-2630/6/1/052 (2004).

Article  ADS  Google Scholar 

Vitale, W. A., Moldovan, C. F., Paone, A., Schüler, A. & Ionescu, A. M. Fabrication of CMOS-compatible abrupt electronic switches based on vanadium dioxide. Microelectron. Eng. 145, 117–119,

doi:10.1016/j.mee.2015.03.055 (2015).

Article  CAS  Google Scholar 

Zimmers, A. et al. Role of Thermal Heating on the Voltage Induced Insulator-Metal Transition in VO2. Phys. Rev. Lett. 110, 56601, doi:10.1103/PhysRevLett.110.056601 (2013).

Article  ADS  CAS  Google Scholar 

Freeman, E. et al. Nanoscale structural evolution of electrically driven insulator to metal transition in vanadium dioxide. Appl. Phys. Lett. 103, 263109, doi:10.1063/1.4858468 (2013).

Article  ADS  Google Scholar 

Hu, C.-L. Self-sustained oscillation in an Rh-C or Rh-L circuit containing a hysteresis resistor Rh. IEEE Trans. Circuits Syst. 33, 636–641, doi:10.1109/TCS.1986.1085968 (1986).

Article  Google Scholar 

Cavalleri, A., Dekorsy, T., Chong, H. H. W., Kieffer, J. C. & Schoenlein, R. W. Evidence for a structurally-driven insulator-to-metal transition in VO2: A view from the ultrafast timescale.

Phys. Rev. B 70, 161102, doi:10.1103/PhysRevB.70.161102 (2004).

Article  ADS  Google Scholar 

Kübler, C. et al. Coherent Structural Dynamics and Electronic Correlations during an Ultrafast Insulator-to-Metal Phase Transition in VO2. Phys. Rev. Lett. 99, 116401,

doi:10.1103/PhysRevLett.99.116401 (2007).

Article  ADS  PubMed  Google Scholar 

Kar, A. et al. Intrinsic electronic switching time in ultrathin epitaxial vanadium dioxide thin film. Appl. Phys. Lett. 102, 72106, doi:10.1063/1.4793537 (2013).

Article  Google Scholar 

Davis, C. The rotation of eigenvectors by a perturbation. J. Math. Anal. Appl. 6, 159–173, doi:10.1016/0022-247X(63)90001-5 (1963).

Article  MathSciNet  MATH  Google Scholar 

Wolpert, D. H. & Macready, W. G. No free lunch theorems for optimization. IEEE Trans. Evol. Comput. 1, 67–82, doi:10.1109/4235.585893 (1997).

Article  Google Scholar 

Culberson, J. & Gent, I. Frozen development in graph coloring. Theor. Comput. Sci. 265, 227–264, doi:10.1016/S0304-3975(01)00164-5 (2001).

Article  MathSciNet  MATH  Google Scholar 

Mammen, D. L. & Hogg, T. A New Look at the Easy-hard-easy Pattern of Combinatorial Search Difficulty. J Artif Int Res 7, 47–66, doi:10.1080/14786419.2014.986658 (1997).

MathSciNet  MATH  Google Scholar 

Vlasie, R. D. Systematic generation of very hard cases for graph 3-colorability. In Seventh International Conference on Tools with Artificial Intelligence, 1995. Proceedings 114–119,

doi:10.1109/TAI.1995.479412 (1995).

Cheeseman, P., Kanefsky, B. & Taylor, W. M. Where the Really Hard Problems Are. In Proceedings of the 12th International Joint Conference on Artificial Intelligence - Volume 1 331–337

(Morgan Kaufmann Publishers Inc., 1991).

Brélaz, D. New Methods to Color the Vertices of a Graph. Commun ACM 22, 251–256, doi:10.1145/359094.359101 (1979).

Article  MathSciNet  MATH  Google Scholar 

Cliques, coloring, and satisfiability: second DIMACS implementation challenge, October 11–13, 1993 (American Mathematical Society, 1996).

Download references

A.P. and M.J. were supported by the National Science Foundation under grant 1640081, and the Nanoelectronics Research Corporation (NERC), a wholly-owned subsidiary of the Semiconductor

Research Corporation (SRC), through Extremely Energy Efficient Collective Electronics (EXCEL), an SRC-NRI Nanoelectronics Research Initiative under Research Task IDs 2698.001 and 2698.002.

N.S., M.J., and S.D. would also like to acknowledge the support of the Office of Naval Research through award N00014-11-1-0665 and Intel Corporation through a customized Semiconductor

Research Corporation project at the University of Notre Dame. This work was also supported, in part, by the Center for Low Energy Systems Technology (LEAST), one of six centers of STARnet, a

Semiconductor Research Corporation program sponsored by MARCO and DARPA.

Author informationAuthors and Affiliations Georgia Institute of Technology, Atlanta, GA, USA

Abhinav Parihar & Arijit Raychowdhury

University of Notre Dame, Notre Dame, IN, USA

Nikhil Shukla, Matthew Jerry & Suman Datta

AuthorsAbhinav PariharView author publications You can also search for this author inPubMed Google Scholar

Nikhil ShuklaView author publications You can also search for this author inPubMed Google Scholar

Matthew JerryView author publications You can also search for this author inPubMed Google Scholar

Suman DattaView author publications You can also search for this author inPubMed Google Scholar

Arijit RaychowdhuryView author publications You can also search for this author inPubMed Google Scholar

A.P. worked on the development of the theory, simulation frameworks and mathematical models. N.S. and M.J. worked on the experiments. A.R. advised A.P. and participated in the problem

formulation. S.D. advised N.S. and M.J. and also participated in the design of experiments and problem formulations.

Corresponding authors Correspondence to Abhinav Parihar or Arijit Raychowdhury.

The authors declare that they have no competing interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

TextRights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or

format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or

other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in

the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the

copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this articleCite this article Parihar, A., Shukla, N., Jerry, M. et al. Vertex coloring of graphs via phase dynamics of coupled oscillatory networks. Sci Rep 7, 911 (2017).

https://doi.org/10.1038/s41598-017-00825-1

Download citation

Received: 04 October 2016

Accepted: 14 March 2017

Published: 19 April 2017

DOI: https://doi.org/10.1038/s41598-017-00825-1

Share this article Anyone you share the following link with will be able to read this content:

Get shareable link Sorry, a shareable link is not currently available for this article.

Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative

This article is cited by A CMOS-compatible oscillation-based VO2 Ising machine solver Olivier MaherManuel JiménezSiegfried Karg Nature Communications (2024)

Mott neurons with dual thermal dynamics for spatiotemporal computing Gwangmin KimJae Hyun InKyung Min Kim Nature Materials (2024)

Computing with oscillators from theoretical underpinnings to applications and demonstrators Aida Todri-SanialCorentin DelacourFilip Sabo npj Unconventional Computing (2024)

Scalable almost-linear dynamical Ising machines Aditya ShuklaMikhail ErementchoukPinaki Mazumder Natural Computing (2024)

High speed universal NAND gate based on weakly coupled RF MEMS resonators Mahdi AttarReza Askari Moghadam Microsystem Technologies (2024)