In order to make the WA problem compatible with the SimCIM quantum-inspired optimization algorithm [35], we first consider a transformation, allowing one to convert the IP problem Eqs 3–5 into a QUBO form as follows: s T Q s → min (6) for a certain binary vector s and the symmetric real matrix Q. This problem is equivalent to finding a configuration of binary-state particles (“spins”) that minimizes the energy: H s = s T Q s , (7) where the Ising Hamiltonian H consists of only single-order terms (energies of individual spins in an external magnetic field) and pair-wise interactions between spins. Although spin variables usually are considered to take values ±1, the transition to a binary form is quite straightforward [13].

A known way [7] to transform a graph coloring problem to the QUBO form is to set s ≔ x (here we treat x as a N _V W-dimensional vector) and use the Hamiltonian of the form: H x = H 1 x + H 2 x , (8) where H 1 x = ∑ v = 1 N V ( 1 − ∑ i = 1 W x v i ) 2 , (9) H 2 x = ∑ u , v ∈ E ∑ i = 1 W x u i x v i . (10)

One can see that H 1 ( x ) > 0 in the case where a single node is assigned with two distinct colors, while H 2 ( x ) > 0 when two adjacent vertices are assigned the same color. If minimization routine provides some x such that H ( x ) = 0 , then x defines a correct coloring with at most W colors. Therefore, an ability to solve the QUBO problem corresponding to the Hamiltonian (Eq. 8) guarantees one to solve a decision problem of whether it is possible to color a graph with at most W colors. Since it is always possible to color a graph with W = N _V colors, a minimal number of colors can be obtained, for example, by using a standard binary search with at most ⌈ log₂(N _V )⌉ iterations. We note that this approach is quite sensitive to possible imperfections of QUBO problem solutions, especially at first iterations of the binary search. An alternative way is to decrease W by a unit at each step that, however, results in a possible increase of iteration numbers up to O ( W start ) , where W _start is the initial upper bound for a color number.

We propose an improved approach for solving a graph coloring problem by developing an alternative transformation into a QUBO form. In our approach, we pursue two major goals. The first is decreasing the number of QUBO problems to be solved. The second is making the whole algorithm robust against the possibility of finding not optimal, but some suboptimal solutions for a particular QUBO problem. We note that these points are of particular importance in the framework of using (quantum-inspired) annealing for solving QUBO problems.

The main idea of our approach is to consider an extended N _V (W+1)-dimensional binary vector s ≔( w , x ) and take the target Hamiltonian in the following form: H w , x = c 0 H 0 w + c 1 H 1 x + H 2 x + c 2 H 3 w , x , (11) where H 0 w = ∑ i = 1 W w i , (12) H 3 w , x = ∑ u , v ∈ E ∑ i = 1 W 1 − w i x u i + x v i , (13) and c _i are positive coefficients satisfying a particular constraint (see more details in Section 5.1). Minimization of this Hamiltonian provides us the solution vector ( w , x ) such that the optimal number of wavelength is encoded in w by non-zero values. We note that the term H 0 ( w ) grows with the total number of used wavelengths; H 1 ( x ) and H 2 ( x ) have the same form as in Eq. 8; H 3 ( w , x ) is responsible for the relationship w _i ≥ x _vi , which becomes positive when the relation is violated. Both terms H 2 ( x ) and H 3 ( w , x ) correspond to inequalities in Eq. 5 in the IP form (see Section 5).

The complete algorithm of solving a graph coloring problem (WA problem) is shown in Algorithm 1. The algorithm uses a subroutine m a k e _ q u b o ( G , W ) that generates the corresponding QUBO matrix Q with respect to the Hamiltonian (Eq. 11), given the input graph G and the target number of the wavelengths W. The QUBO problem is then solved with the subroutine s o l v e _ q u b o ( Q ) , which finds the optimal spin vector s = ( w , x ) using the quantum-inspired SimCIM approach for the QUBO matrix Q, as defined in [35]. In order to check the validness of the obtained solution, we use c h e c k _ c o l o r i n g ( G , x ) that validates the fulfillment of Eqs 4, 5.

Algorithm 1

Solving graph coloring problem with improved transformation:

One can see that, if s o l v e _ q u b o ( Q ) provides an optimal solution, then the whole problem is solved in the first iteration. However, even in the case when the obtained solution is suboptimal, the updated problem with the reduced upper bound W becomes easier to solve, and the algorithm converges with a few numbers of iterations.

In this study, we solve the WA problem and obtain results with the use of 1) the proposed technique based on quantum-inspired optimization SimCIM [35] (with the improved approach, see Section 5), 2) industry-grade commercial Gurobi optimization software, and 3) the open-source mixed-integer programming solver—GLPK. We note that in the case of the quantum-inspired optimization with SimCIM, we solve the problem in the QUBO form (Eq. 11), whereas in the case of Gurobi and GLPK, we use the IP formulation of graph coloring [see (Eqs 3–5)]. Additionally, we include the results obtained via the largest-degree-first (LDF) heuristics used as the baseline since it allows one to instantly produce feasible coloring without numerical optimization. We also ran the experiments for original QUBO transformation proposed in [7] and compared them to our proposed QUBO in the Table A1.

Our numerical experiments have been performed on a synthetic dataset of 900 randomly generated graphs with varying node numbers and edge probabilities (for details, see Section 5.3). The main characteristics that we are interested in are time-to-solution (TTS) and the number of colors in the obtained solution. The total run time has been limited by 300 s, and the best solutions have been compared. Results are averaged over 90 runs for each graph size (for details, see Table 1). For all numerical experiments, we use the same hardware set, which is based on Xeon E-2288G 3.7 GHz CPU, 128GB RAM, and GeForce GTX1080 8 GB graphics card.

Our results indicate that the quantum-inspired technique SimCIM demonstrates a behavior comparable with Gurobi in the case of small nodes (10–30 nodes). Moreover, the run time of SimCIM is better for large-scale (90 and 100 nodes) graphs, as it is indicated in Table 2. Such a trend can be explained. As the number of nodes increases, the number of inequalities in the ILP formulation of the problem grows rapidly. The number of inequalities is equal to the product of the number of edges by the number of colors available for coloring the vertices of the graph. So, the complexity of the problem for the ILP solver increases rapidly with the number of nodes. GLPK shows a stable result up to 30 nodes and becomes unstable further after a timeout interruption without any solution with more than 10 percent instances. We note that the comparison between our quantum-inspired approach and Gurobi is conducted in the common setting, so its additional tuning for obtaining better results is also possible. At the same time, we find it interesting that the quantum-inspired technique shows comparable or superior results in harder, industry-relevant, combinatorial optimization problem.

Although our goal was to demonstrate the applicability of a quantum-inspired graph coloring algorithm for the wavelength assignment problem, our approach can be applied to a variety of problems, in particular from the field of scheduling [42].

Assuming that we have the set of jobs to schedule, every job requires one time slot and some jobs cannot be executed at the same time due to some interference with each other; we need to determine the minimal time when every job will be finished or how many time slots they will occupy. One can build the graph so that vertices correspond to the jobs, and two vertices are connected if these jobs cannot be executed at the same time. The colors of vertices represent the time slots to assign; so, a graph has k number of colors if the jobs can be executed in k time slots.

Using our approach, we take the proposed Hamiltonian in Eq. 11 and redefine its variables so that the following expression is obtained: x v i = 1 , if vertex  v  is assigned time slot  i , 0 , otherwise , (14) and w i = 1 , if  i - th time slot is assigned , 0 , otherwise . (15)

That way, the job-scheduling problem can be solved using quantum-inspired annealing analogously to the WA problem.

The same approach can be implemented for tasks from other fields, such as computer register allocation [43], storage of chemicals [44], and printed circuit board testing [45].

The main step in solving an optimization problem using the quantum and quantum-inspired annealing method is to map the problem of interest to the energy Hamiltonian (so-called Ising Hamiltonian) so that the quantum device could find the ground state that corresponds to the optimum value of the objective function. In this study, we formulate a mapping of the graph coloring problem into the QUBO form given by Eq. 6. There is a well-known transformation of the graph G = (V, E) coloring decision model [7] that shows the possibility of coloring with a constant number of colors W, but we represent novel QUBO transformation that could minimize the number of colors and implement the original problem statement (Eqs 3–5).

The objective function ∑ i = 1 W w i could be exactly mapped to the QUBO form: H 0 w = ∑ i = 1 W w i , (16) where w = (w ₁, ..., w _W ) is a binary vector indicating colors used in coloring. The constraint ∑ i = 1 W x v i = 1 for every v ∈ V after mapping takes the form H 1 x = ∑ v = 1 N V 1 − ∑ i = 1 W x v i 2 , (17) where N _V is the number of nodes in G.

The situation with the second constraint x _ui + x _vi ≤ w _i for every i ∈ {1, ..., W} and (u, v) ∈ E appears to be more complicated. One can see that it involves three variables and thus cannot be directly embedded into a two-body Hamiltonian. However, we can use the following trick. One can easily check that for arbitrary a, b, c ∈ {0, 1}, the following equivalence holds: a + b ≤ c ⇔ a b = 0 , 1 − c a + b = 0 . (18) This fact allows us to embed the conditions x _ui + x _vi ≤ w _i into two Hamiltonians: H 2 x = ∑ u , v ∈ E ∑ i = 1 W x u i x v i , (19) H 3 w , x = ∑ u , v ∈ E ∑ i = 1 W 1 − w i x u i + x v i . (20)

The resulting Hamiltonian consists of all components’ sum: H w , x = c 0 H 0 w + c 1 H 1 x + H 2 x + c 2 H 3 w , x , (21) where c ₀, c ₁, and c ₂ are positive constants standing for a positive penalty value. We note that the sum H 1 ( x ) + H 2 ( x ) is exactly matched with the classical decision problem [7] and responsible for the correct coloring of the graph. Therefore, H 1 ( x ) and H 2 ( x ) are grouped with the same penalty coefficient c ₁. Coefficients c ₀, c ₁, and c ₂ should be set manually using the following criteria: the penalty value c ₁ should be high enough to the keep the final solution from violating constraints. At the same time, a very big penalty value can overwhelm the target function, making it difficult to distinguish solutions of different qualities. We establish inequalities for constraint coefficients that show the equivalence of IP and QUBO models of a problem.

In this study, we demonstrate how to construct an operator matrix Q of our QUBO model for the WA problem. Recall that we take the binary vector of the QUBO problem in the form s = ( w , x ), i.e., enumerate K = (N _V +1)W and binary variables s _k and link them to our model variables as follows: s k = w k , k = 1 , … , W , x u i , k = u W + i , (34) where u = 1, … , N _V , i = 1, … , W.

The goal is to find the vector s that minimizes the quadratic form s ^T Q s , and we show that it is equivalent to minimizing energy of the Hamiltonian (Eq. 11). Let us denote A the adjacency matrix of the network graph G = (V, E) so that a _uv = 1 if (u, v) ∈ E, and a _uv = 0 otherwise. We note that the sum of the vth column of A equals the degree of the vertex v, and the sum of all vertex degrees is 2N _E . We rewrite the operator (Eq. 11) terms H 0 ( w ) , H 1 ( x ) , H 2 ( x ) , and H 3 ( w , x ) as follows: H 0 w = ∑ i = 1 W w i 2 , (35) H 1 x = ∑ v = 1 N V 1 − ∑ i = 1 W x v i 2 = ∑ v = 1 N V ∑ i = 1 W x v i 2 − 2 ∑ i = 1 W x v i + N V = ∑ v = 1 N V ∑ i , j = 1 W x v i x v j − 2 ∑ i = 1 W x v i 2 + N V , (36) H 2 x = ∑ u , v ∈ E ∑ i = 1 W x u i x v i = ∑ u , v = 1 N V ∑ i = 1 W a u v x u i x v i , (37) H 3 w , x = ∑ u , v ∈ E ∑ i = 1 W 1 − w i x u i + x v i = ∑ i = 1 W 1 − w i ∑ v = 1 N V d v x v i = − ∑ v = 1 N V ∑ i = 1 W d v w i x v i + ∑ v = 1 N V d v ∑ i = 1 W x v i . (38) In expanding the expression for H 1 ( x ) , we exploit the fact that since x _vi is binary, then x v i 2 = x v i . Also, we note that if H 1 ( x ) = 0 , then the last term in H 3 ( w , x ) equals 2N _E .

Considering the equalities (Eqs 35–38) for Hamiltonian terms H 0 ( x ) , H 1 ( x ) , H 2 ( x ) and H 3 ( w , x ) , we construct a QUBO operator as a block matrix as follows: Q = Q 11 Q 12 Q 21 Q 22 , (39) where Q 11 = c 0 E W , (40) Q 12 = − c 2 2 D ⊗ E W , Q 21 = Q 12 T , (41) Q 22 = c 1 E N V ⊗ I W − 2 E W + c 1 A ⊗ E W . (42)

Here, E _W denotes the identity matrix of size W, I _W denotes a matrix with all elements equal to 1 of those of size W, and D = ( d 1 , . . . , d N V ) is a row vector of graph vertex degrees. We also employ the fact that the terms of the form ∑ u , v = 1 N V ∑ i , j = 1 W c u v h i j x u i x v j , (43) for some coefficients c _uv = c _vu and h _ij = h _ji can be represented by a quadratic form defined by the Kronecker product C ⊗ H, where C and H are matrices of c _uv and h _ij , respectively. Matrix Q is constructed so that the Q ₁₁ submatrix corresponds to the term H 0 ( x ) of the Hamiltonian (Eq. 11), the Q ₁₂ submatrix is for H 3 ( w , x ) , and Q ₂₂ is for H 1 ( x ) + H 2 ( x ) .

It is worth emphasizing that it is the structure of the encoding problem parameters into the spin vector, given by Eq. 34, that allows us to represent submatrices Q ₁₂, Q ₂₁, and Q ₂₂ in the form of Kronecker products. This feature of QUBO submatrices significantly speeds up their assembly using standard mathematical packages, e.g., n u m p y and s c i p y .

We generate datasets that are used in binomial graphs [47], or Erdös–Rényi graphs, which have two parameters for generation: the number of nodes N _V and the probability of an edge occurrence p. Each of possible N = N _V ⋅ (N _V −1)/2 edges is chosen with probability p. The number of edges N _E is drawn randomly from the binomial distribution: P N E = x =   x N p x ⋅ q N − x . (44)

To take into account sparse and dense graphs, various probability p options from .1 to .9 with an interval of .1 have been chosen; the number of graph nodes has been varied 10 to 100 with a step of 10. For each pair (n, p), 10 connected graphs have been generated with different seed parameters. We note that disconnected graphs are not included the dataset. The overall characteristics of the dataset are given in Table 3.

Optimal penalty values guarantee the fulfillment of constraints for an optimal solution, but large values of c ₁ and c ₂ reduce the contribution of the initial objective function to the total energy and significantly increase the time to find the optimal solution. Our approach to solve this problem is as follows:

1. Set the minimum possible penalty values c ₁ and c ₂ using trial runs so that the contribution of the objective function is sufficient.

The following penalty values were set for the tests: c 0 = 1 , c 1 = 10 + p N V , c 2 = 2.5 . (45)

SimCIM [35] is an example of a quantum-inspired annealing algorithm, which works in an iterative manner. SimCIM can be used for sampling low-energy spin configurations in the classical Ising model, and the Hamiltonian can be written as follows: H = ∑ i h i s i + ∑ < i , j > J i j s i s j , (46) where J represents the spin–spin interaction, h represents the external field, and s _i are the individual spins on each of the lattice sites. The Ising Hamiltonian can be directly transformed to a QUBO problem [13], and then, quantum annealing can be applied to any optimization problem, which can be expressed into the Quadratic Unconstrained Binary Optimization (QUBO) form. The SimCIM algorithm treats each spin value as a continuous variable s _i ∈ [−1, 1]. Each iteration of the algorithm starts with calculating the mean field of the following form: Φ i = ∑ j ≠ i J i j s j + h i , (47) which acts on each spin by all other spins. Then, the gradients for the spin values are calculated as follows: Δ s i = p t s i + ζ Φ i + N 0 , σ , (48) where p _t is a dynamic parameter dependent on the SimCIM annealing process, the overall feed forward factor is ζ, and N (0, σ) is a random variable sampled from the Gaussian distribution with zero mean and standard deviation σ. Then, the spin values are updated according to s _i ← ϕ(s _i +Δs _i ), where ϕ(x) is the activation function: ϕ x = x   for  | x | ≤ 1 ; x / | x | , o t h e r w i s e . (49)

After multiple updates, the spins will tend to either −1 or +1, and the final discrete spin configuration is obtained by taking the sign of each s _i .

In our implementation, we added several improvements to the SimCIM algorithm defined in the original paper [35]. In particular, we normalized the value of the Gaussian noise to a gradient norm and introduced gradient quantization, which made the solver more stable near optimum points.