Quantum-inspired optimization for wavelength assignment

Problems related to wavelength assignment (WA) in optical communications networks involve allocating transmission wavelengths for known transmission paths between nodes that minimize a certain objective function, for example, the total number of wavelengths. Playing a central role in modern telecommunications, this problem belongs to NP-complete class for a general case, so that obtaining optimal solutions for industry relevant cases is exponentially hard. In this work, we propose and develop a quantum-inspired algorithm for solving the wavelength assignment problem. We propose an advanced embedding procedure for this problem into the quadratic unconstrained binary optimization (QUBO) form having an improvement in the number of iterations with price-to-pay being a slight increase in the number of variables ("spins"). Then we compare a quantum-inspired technique for solving the corresponding QUBO form against classical heuristic and industrial combinatorial solvers. The obtained numerical results indicate on an advantage of the quantum-inspired approach in a substantial number of test cases against the industrial combinatorial solver that works in the standard setting. Our results pave the way to the use of quantum-inspired algorithms for practical problems in telecommunications and open a perspective for the further analysis of the employ of quantum computing devices.


I. INTRODUCTION
Optimization is a tool with applications across various technologies [1]. However, solving complex real-world optimization problems is computationally intensive even in the case of using advanced, specialzed hardware. Quantum computers are widely believed to be useful for solving computationally difficult optimization problems beyond the capability of existing computing devices is to use quantum optimization [2][3][4][5][6]. A general approach consists in encoding a cost function in a quantum Hamiltonian [7], so that its low-energy state is obtained starting from a generic initial state. Among existing methods to achieve such dynamics, quantum annealing offers physical implementations of a non-trivial size [8]. Quantum annealing is by now explored for analysis of various areas, such as chemistry calculations [9,10], lattice protein folding [11,12], genome assembly [13,14], solving polynomial systems of equations for engineering applications [15] and linear equations for regression [15], portfolio optimization [16][17][18][19], forecasting crashes [20], finding optimal trading trajectories [21], optimal arbitrage opportunities [22], optimal feature selection in credit scoring [23], foreign exchange reserves management [24], traffic optimization [25][26][27], scheduling [28][29][30][31][32][33], railway conflict management [32,33], and many others [5]. Advances also include the recent experimental demonstration of a superlinear quantum speedup in finding exact solutions for the hardest maximum independent set graphs [34].
Although quantum optimization algorithms suggest an intriguing possibility to solve computationally difficult problems beyond the capability of classical computers, exiting conceptual and technical limitations make it challenging to use it for solving problems of industry relevant sizes. Attempts to simulate quantum computations classically resulted in a new class of algorithms and techniques know as quantum-inspired [35,36].
As soon as these algorithms are compatible with currently existing (classical) hardware, analyzing their limiting capabilities and advantages over classical approaches are required towards their use in practice. Specifically, a way to solve combinatorial optimization problems via simulating the coherent Ising machine (SimCIM) has been proposed [35]. SimCIM algorithm is able to solve optimization problems that are formulated in the quadratic unconstrained binary optimization (QUBO) / Ising form, which can be done for various practically relevant cases [7]. The SimCIM approach has demonstrated capabilities to outperform bona fide coherent Ising machine and existing classical methods for certain GSet graphs. However, one of the arising questions is related to the need in to tune hyperaparameters [35]. For a wide range of benchmark of quantum-inspired heuristic solvers for quadratic unconstrained binary optimization, namely D-Wave Hybrid Solver Service, Toshiba Simulated Bifurcation Machine, Fujitsu Digital Annealer, and simulated annealing on a personal computer, see also Ref. [37].
Design of optical communication network is a specific industrial avenue, in which combinatorial optimisation in ubiquitous. Examples of tasks include finding optimal transmission and reservation paths, frequency allocation, network throughput maximization and many others [38,39]. A notable example is the routing and wavelength assignment (RWA) problem, which consists in allocating transmission wavelengths and finding transmission paths between nodes that minimize the total number of wavelengths. Conventional techniques, such as linear programming and mixed integer programming, are useful for most of the cases; however, the combinatorial nature and hardness of the problems make it extremely challenging to apply these techniques for large-scale problems. It is then reasonable to assume that telecommunication industry may benefit from the use of quantum-inspired algorithm in the near-term horizon and quantum com- FIG. 1. Illustration of the approach. A linear network with generated requests and paths consisting of 5 nodes, 4 edges, and 5 traffic paths is considered: Solid lines represent original edges, and the arrows lines represent traffic paths. One can reduce the WA problem to a graph coloring problem with a simple graph transformation (bottom of the figure): Each traffic path is now considered a vertex; if two traffic paths share (at least) one fiber they are connected by an edge.
In this study, we consider the variant of the RWA problem. To explain more precisely, we focus on the wavelength assignment task for known routes which we further refer to as the wavelength assignment (WA) problem. This problem is generally NP-hard, so its solution is computationally challenging for large sizes. We propose an original way to transform the WA problem to the QUBO form, which makes it compatible with quantum-inspired optimization algorithm and, in principle, quantum annealing hardware. For solving this problem, we develop a technique based on the SimCIM quantum-inspired optimization solver [35] with the use of the Lagrange multipliers for minimizing the number of hyperparameters. Our numerical results indicate on an advantage of the quantum-inspired solver in a number of test cases against the industrial combinatorial solver working on the standard settings.

II. WAVELENGTH ASSIGNMENT PROBLEM (WA)
Let us consider a network connecting a number of endpoints with optical links (see an example in Fig. 1 receiver. A single optical link can be shared between several paths given that each path is assigned different wavelengths. Each path is indicated by the path ID, which uniquely identifies a pair of transmitting/receiving nodes, sequence of interconnecting nodes, and the wavelength ID. The WA problem implies allocation of the wavelength IDs for paths that are pre-computed and known a priori in such a way to meet the target objective, for example, the number of the used wavelengths is minimized 1 . Formally, WA is considered to be correct if and only if it satisfies the following requirements: (i) each path should use a single wavelength and (ii) several paths sharing the same edge should have different wavelengths.
The problem of finding correct wavelength allocation under given constraints is equivalent to the coloring problem [7] in transformed graph G = (V, E), where nodes V and edges E representing paths and their intersections in fibers, correspondingly (two nodes from V are connected if and only if the corresponding paths have an intersection within the optical network). Let N V and N E denote numbers of vertices and edges of G, respectively. Later we interchangeably use terms wavelengths and colors since the underlying problems are formally identical. The example of the correspondence of network paths to graph coloring mapping is shown in Fig. 2.
In order to define a particular coloring of graph G with at most W colors, we introduce a two collections of aux-iliary variables. The first variable is x that consists of N V W binary variables The second one, denoted by w, consists of W binary variables Employing x and w, the problem of finding an correct allocation with minimum number of the used wavelengths not exceeding some maximal number W ≥ 1, can be formulated as an integer programming (IP) problem of the following form: One can see that constraint (4) assures that each vertex is assigned to exactly one wavelength, while constraint (5) indicates that two adjacent vertices are not assigned the same wavelength. This problem is generally NP-hard, so its solution is computationally challenging for large sizes. As it is shown below, the QUBO reduction makes the problem compatible with quantum-inspired algorithms that can shift tractability boundaries to higher problem sizes. While such reduction usually involves additional overheads in the problem size due to auxiliary variables, the overheads can be compensated by the computational advantage of quantum-inspired solvers leading to better overall results.

A. Transforming the WA problem to a QUBO form
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 (3)-(5) into a QUBO form as follows: 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 where the Ising Hamiltonian H consists of only singleorder terms (energies of individual spins in 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 where One can see that H 1 (x) > 0 in the case where 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 Hamiltonian (8) garantees 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 2 (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 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 colors number.

B. Improving QUBO transformation for quantum-inspired annealing
We propose an improved approach for solving 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 solution for a particular QUBO problem. We note that these points are of particular importance in the framework of using (quantuminspired) 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: and c i are positive coefficients satisfying a particular constraint (see more details in Methods V A). Minimization of this Hamitonian 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); and 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 (5) in the IP form (see Methods V). The complete algorithm of solving graph coloring problem (WA problem) is shown in Algorithm 1. The algorithm employs a subroutine make qubo(G, W ) that generates the corresponding QUBO matrix Q with respect to Hamiltonian (11), given input graph G and the target number of the wavelengths W . The QUBO problem is then solved with subroutine solve qubo(Q), which finds the optimal spins vector s = (w, x) using the quantuminspired SimCIM approach for the QUBO matrix Q as defined in Ref. [35]. In order check the validness of obtained solution, we use check coloring(G, x) that validates the fulfilment of Eq. (4) and Eq. (5).

Algorithm 1 Solving graph coloring problem with improved transformation
Require: W is the initial upper bound on the number of wavelengths One can see that, if solve qubo(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 sub-optimal, the updated problem with the reduced upper bound W becomes easier to solve, and the algorithm converges with a few numbers of iterations.

C. Numerical results
Here we solve the WA problem and obtain results with use of (i) the proposed technique based on quantuminspired optimization SimCIM [35] (with the improved approach, see Methods), (ii) industry grade commercial Gurobi optimization software, and (iii) 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 (11), whereas in case of Gurobi and GLPK we use the IP formulation of graph coloring [see Eqs.
(3)- (5)]. Additionally, we include the results obtained via 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 transormation proposed in Ref.
[7] and compared them to our proposed QUBO in the Table IV in Appendix.  Our numerical experiments have been performed on a synthetic dataset of 900 randomly generated graphs with varying nodes number and edge probability (for details, see Methods V C). The main characteristics that we are interested in are time-to-solution (TTS) and the number of colors in the obtained solution. The total runtime has been limited by 300 seconds, and the the best solutions have been compared. Results are averaged over 90 runs for each graph size (for details, see Table V). For all numerical experiments, we use the same hardware set, which is based on Xeon E-2288G 3,7GHz CPU, 128GB RAM, and GeForce GTX1080 8GB graphics card.  Our results indicate that the quantum-inspired technique SimCIM demonstrates behaviour comparable with Gurobi in the case of small (10-30 nodes). Moreover, the runtime of SimCIM is better for large-scale (90 and 100 nodes) graphs as it is indicated in Table II. Such 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 interrupt 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 quantum-inspired technique shows comparable or superior results in harder, industry relevant combinatorial optimization problem.

D. Other potential applications
While our goal was to demonstrate the applicability of quantum-inspired graph coloring algorithm for wavelength assignment problem, our approach can be applied to a variety of problems, in particular from the field of scheduling [42].
Assuming we have the set of jobs to schedule, every job requires one time slot and some jobs can not 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 can't be executed at the same time. The colors of vertices represent time slots to assign, so graph has k number of colors if the jobs can be executed in k time slots.
Using our approach we take proposed Hamiltonian in Eq. (11) and redefine its variables so that and That way the jobs scheduling problem can be solved using quantum-inspired annealing analogously to WA problem.
The same

IV. CONCLUSION
A search for new approaches to solving practicallyrelevant optimization problems is a clear goal for many industrial applications since even minor improvement on a large scale may generate serious economical impact. In this domain, much attention is paid to quantum computing, which is believed to be useful for such class of problems. At the current technological level, practical quantum advantage, for example, in optimization is still needed to be achieved. An interesting part of this research is understanding of the physical origin of the potential advantages of quantum computing technologies. Attempts to simulate quantum computation classically resulted in a new class of algorithms and methods know as quantum-inspired, which are ready to be tested for industry-relevant problems.
In this work, we have considered the industry important problem in the field of telecommunications. We have demonstrated a way to make it compatible with quantum and quantum-inspired techniques. Interestingly, our numerical results have indicated on an advantage of the quantum-inspired solver in a number of test cases against the industrial combinatorial solver working on the standard settings.
One may expect that the additional tuning of the industry grade commercial optimization solver may result in a substantial improving of its performance. At the same time, studying the origins of the advantages of the quantum-inspired approach, which are largely beyond the scope of the present proof-of-concept demonstration, would allow its further progress as well.
We would like to note that our comparison is limited by the upper bound of 100 nodes, since it allows us to run all solvers in equivalent hardware setup using CPU mode on single core. Further analysis of larger graphs requires running SimCIM solver on GPU card, which gives significant acceleration factor not directly available in conventional MIP algorithms, which are heavily dependent on graph processing routines. As for the multi-core CPU execution environment, some MIP solvers can benefit from such setup by running various optimization strategies and hyperparameters simultaneously. Such speed up quickly reaches saturation point at the level of 8 16 cores (with around 2x improvement in accordance with Gurobi experiments, see slide 26 [46]) and demonstrates no substantial improvement at higher concurrency levels. On the other side, quantum-inspired approach exploits parallelism on the level of starting optimization points, which demonstrates slower, but stable performance increase at the higher levels of concurrency (100 ∼ 1000 parallel units of execution). Thus, we conduct our benchmarks exclusively using CPU mode on single-core to avoid bias towards either solution approach. In order to maintain fairness of comparison for larger graphs our benchmark routine should be further revised to account for heterogeneous (CPU/CPU multi-core vs GPU/multi-GPU) computing environments.

A. Hamiltonian of wavelength assignment problem
The main step in solving an optimization problem using quantum and quantum-inspired annealing is to map the problem of interest to the energy Hamiltonian (socalled Ising Hamiltonian), so the quantum device could find the ground state that corresponds to the optimum value of the objective function. Here we formulate a mapping of the graph coloring problem into QUBO form given by Eq. (6). There is a well-known transformation of the graph G = (V, E) coloring decision model [7] that shows possibility of coloring with some constant number of colors W , but we represent novel QUBO transformation that could minimize number of colors and implement original problem statement Eq. (3-5).
The objective function W i=1 w i could be exactly mapped to the QUBO form: where, recall, w = (w 1 , . . . , w W ) is a binary vector indicating colors used in coloring. The constraint W i=1 x vi = 1 for every v ∈ V after mapping takes the form 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 can not 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: This fact allows us to embed the conditions x ui +x vi ≤ w i into two Hamiltonians: The resulting Hamiltonian consist of all components sum: where c 0 , c 1 , and c 2 are positive constants stand 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), H 2 (x) are grouped with the same penalty coefficient c 1 . Coefficients c 0 , c 1 , and c 2 should be set manually, using the following criteria: the penalty value c 1 should be high enough to keep the final solution from violating constraints. At the same time, too 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. Proposition (QUBO penalty coefficients selection). Consider an IP problem given by Eq. (3)-(5) for some maximal colors number W and some graph G = (V, E) with N E edges. If the IP problem has a solution, then the corresponding QUBO problem, given by Hamiltonian (21) with penalty coefficients satisfying has a solution, equivalent to the solution of the IP problem.
Proof. First, let us rewrite Hamiltonian (21) in the form where Note that A, B, and C can take non negative integer values only. Let (w I , x I ) and (w Q , x Q ) be solutions of the IP and QUBO problems correspondingly. Our goal is to prove that (i) i.e., (x Q , w Q ) defines a correct coloring, and (ii) i.e., the solution of the QUBO problem coincides with the one of the IP problem. First, let us see that Eq. (22) assures B(x Q ) = 0. The proof of this part is by a contradiction. Suppose that B(x Q ) ≥ 1. Consider the difference of energy functions . (28) The correctness of the IP solution implies B(x I ) = 0, and so B(x Q ) − B(x I ) ≥ 1. The differences in terms with A and C can be lower bounded by the corresponding extreme values: In this way, Eq. (28) transforms into given constraint (22). However, this result contradicts with the fact that (w Q , x Q ) provides the minimal energy. Therefore, B(x Q ) = 0, and We then prove that provided C(w I , x I ) = 0 and the second constraint (23).
Finally, A(w Q ) = A(w I ), since otherwise, either there exist a solution for the QUBO problem that is better than (w Q , x Q ), or (x I , w I ) is not the true solution the IP problem.
Therefore, the optimal solution to the QUBO problem appears to be equivalent to the optimal solution to the corresponding IP problem.

B. Wavelength assignment QUBO transformation
Here 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 binary variables s k and link them to our model variables as follows: where u = 1, . . . , N V , i = 1, . . . , W .
The goal is to find vector s that minimizes quadratic form s T Qs and we show that it is equivalent to minimizing energy of Hamiltonian (11). Let us denote A the adjacency matrix of network graph G = (V, E) so that a uv = 1 if (u, v) ∈ E and a uv = 0 otherwise. We note that sum of v-th column of A equals the degree of the vertex v and the sum of all vertex degrees is 2N E . We rewrite operator (11) terms H 0 (w), H 1 (x), H 2 (x) and H 3 (w, x) as follows: In expanding the expression for H 1 (x) we exploit the fact that since x vi is binary then x 2 vi = x vi . Also we note that if H 1 (x) = 0 then the last term in H 3 (w, x) equals 2N E .
Considering the equalities (35-38) for Hamiltonian terms H 0 (x), H 1 (x), H 2 (x) and H 3 (w, x), we construct QUBO operator as block matrix as follows: where 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 terms of the form for some coefficients c uv = c vu and h ij = h ji can be represented by a quadratic form defined by Kronecker product C⊗H, where C and H are matrices of c uv and h ij correspondingly. One can that matrix Q is constructed so that Q 11 submatrix corresponds to the term H 0 (x) of Hamiltonian (11), Q 12 submatrix is for H 3 (w, x) and It is worth to emphasize that it is the structure of encoding problem parameters into the spin vector, given by (34), that allow us to represent submatrices Q 12 , Q 21 , and Q 22 in the form of Kronecker products. This feature of QUBO submatrices significantly speeds up their assembly using standard mathematical packages, e.g., numpy and scipy.

C. Dataset generation
We generate dataset are used 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. Number of edges N E are drawn randomly from binomial distribution: To take into account sparse and dense graphs, various probability p options from 0.1 to 0.9 with an interval of 0.1 have been chosen, the number of graph nodes has been varied from 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 the Table III

D. Setting penalty values
Optimal penalty values guarantee the fulfilment of constraints for an optimal solution, but large values of c 1 and c 2 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 1 and c 2 using trial runs, so that the contribution of the objective function is sufficient; 2. Use all SimCIM iterations to select feasible solutions; 3. Take the feasible solution with the lowest energy.
E. Quantum-inspired annealing using SimCIM 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 which Hamiltonian can be written as: where J represents the spin-spin interaction, h represents the external field, and the 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 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: which act on each spin by all other spins. Then the gradients for the spin values are calculated as follows: where p t is a dynamic parameter dependent on Sim-CIM annealing process, overall feedforward factor ζ 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 x for |x| ≤ 1; x/|x|, otherwise.
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 add several improvements to SimCIM algorithm defined in the original paper [35]. In particular, we normalize the value of the Gaussian noise to gradient norm and introduced gradient quantization, which made the solver more stable near optimum points. Average results (lower is better)