Scalable Node-Disjoint and Edge-Disjoint Multi-wavelength Routing

Probabilistic message-passing algorithms are developed for routing transmissions in multi-wavelength optical communication networks, under node and edge-disjoint routing constraints and for various objective functions. Global routing optimization is a hard computational task on its own but is made much more difficult under the node/edge-disjoint constraints and in the presence of multiple wavelengths, a problem which dominates routing efficiency in real optical communication networks that carry most of the world's Internet traffic. The scalable principled method we have developed is exact on trees but provides good approximate solutions on locally tree-like graphs. It accommodates a variety of objective functions that correspond to low latency, load balancing and consolidation of routes, and can be easily extended to include heterogeneous signal-to-noise values on edges and a restriction on the available wavelengths per edge. It can be used for routing and managing transmissions on existing topologies as well as for designing and modifying optical communication networks. Additionally, it provides the tool for settling an open and much debated question on the merit of wavelength-switching nodes and the added capabilities they provide. The methods have been tested on generated networks such as random-regular, Erd\H{o}s R\'{e}nyi and power-law graphs, as well as on the UK and US optical communication networks. They show excellent performance with respect to existing methodology on small networks and have been scaled up to network sizes that are beyond the reach of most existing algorithms.


I. INTRODUCTION
Optical communication networks underpin the global digital communications infrastructure and carry most of the Internet traffic. They comprise thousands of kilometers of optical fibers, organized in a complex web of constituent sub-networks including the Internet backbone, Metro access and Data Center farms. The exponential growth in Internet traffic and energy consumption threatens to overload the existing infrastructure and a capacity crunch is looming [1]. Not only that a matching growth in infrastructure is infeasible, it raises fundamental questions on the ultimate capacity of optical communication networks and the manner in which we could optimize their use. The next-generation digital infrastructure has to offer flexibility, low latency, high network throughput and resilience.
One of the key requirements is the routing and wavelength assignment (RWA) for all traffic demands across this complex heterogeneous network in a way that optimizes a given objective function, be it low latency, high throughput or resilience [2]. Each optical fiber carries information using light of one or many wavelengths. The latter uses, among others, dense wavelength-division multiplexing (DWDM) methods that employ as many as 80-160 channels of different laser wavelengths [3]. The main constraint in the RWA is that any complete individual route, from source to destination, uses the same single wavelength and that two separate routes using the same wavelength cannot shares the same fiber. This constraint makes the corresponding mathematical problem hard to solve in general.
Route optimization in optical communication networks can be mapped onto the hard computational problem of edge-disjoint routing on a graph, where transceivers (transmitter-receiver) are mapped to vertices (or nodes) and fibers to edges (or links). Given that routes are constrained to be contiguous and interaction between paths is non-localized, local optimization methods are insufficient and global optimization is required. Globally optimal routing of multiple messages or vehicles given a general objective function is a computationallyhard constraint satisfaction problems on its own and has been addressed in the physics literature using scalable and distributed message passing approximation techniques, inspired by statistical physics methodology [4][5][6][7]. Moreover, similar techniques have been suggested also for addressing the single-wavelength Node-Disjoint Paths (NDP) [8] and Edge-Disjoint Paths (EDP) [9] problems where multiple paths of different origin-destination pairs cannot share nodes or edges on a graph, respectively. Generally, both optimization tasks are within the class of NP-hard combinatorial problems [10][11][12][13] and the approximation offered by message passing techniques work well. However, the existing methods developed for singlewavelength routing become intractable in the presence of multiple wavelengths, making them inapplicable for realistic scenarios.
It is worthwhile noting that in some extreme cases these hard computational problems become polynomial in the single wavelength case as discussed in [9]. For instance, when the number of origin-destination pairs is low with respect to the systems size [14] and where all origindestination pairs are identical [15,16]. Nevertheless, the general problems of NDP and EDP routing on graphs are computationally hard even in the single wavelength case and a variety of methods have been used to address them in the context of optical communication networks and more general problems. Alongside established methods, such as integer/linear programming and its variants [17][18][19][20][21] bin-packing based approaches [22], Monte Carlo search [23], post-optimization [24] and greedy algorithms [25][26][27][28], a large number of heuristics have also been used to obtain approximate solutions in both problems, among them genetic algorithms [29][30][31], ant colony optimization [32] and particle swarm optimization [33]. Specifically in the area of optical communication networks, it is common to use integer/linear programming and its variants for small networks, to obtain exact results, and a variety of heuristics for larger networks. In practice, current optical networks use overprovision of capacity to compensate for sub-optimal routing, resulting in both over-engineering and underutilized capacity [34,35].
The main challenge we address here is the RWA under heavy traffic using multiple wavelengths and a very large number of origin-destination pairs under the NDP and EDP constraints and for various objective functions. Globally optimal routing is a non-localized difficult problem but adding the NDP/EDP restrictions and having a large number of different wavelengths increases the complexity considerably, making the problem intractable for existing algorithms [8,9]. We map the globally-optimal routing problem in the presence of multiple wavelengths onto a multi-layer replica of the original graph and utilize probabilistic optimization approaches. The methods developed here are based on message passing techniques, developed independently in several fields including statistical physics, computer science and information theory [36][37][38] but are closely interlinked [39,40]. The methods we develop allow for messages, in the form of conditional probability values to be passed between nodes and the replicated networks representing the different wavelengths, in a way that keeps the algorithms scalable and applicable even for a large number of wavelengths, transmissions (corresponding to source-destination pairs) and nodes.
The main result of this paper is the derivation of principled scalable algorithms, capable of obtaining approximate solutions for routing problems in large graphs, where the number of transmissions is of similar order to that of the number of free variables (vertices/edges) and a large number of wavelengths, under the NDP and EDP constraints and for various objective functions, both convex and concave. The algorithm also accommodate cases where the number of transmissions is much larger than the number of vertices (quadratic with respect to the number of vertices). The computational complexity of the NDP/EDP algorithms for sparse graphs is O(M Q(M + N + Q))/O(M Q(M/N + N + Q)), with N the number of vertices, M the number of transmissions and Q the number of wavelengths. The algorithm has been tested for a variety of sparse network topologies, both synthetic random graphs and real optical commu-nication networks, and for different objective functions, showing excellent results in obtaining high quality approximate solutions. Among the generic networks examined are random regular graphs, Erdős Rényi (ER) [41] and scale free networks [42], while the realistic networks considered include the British 22 nodes (BT22) and US 60 node (CONUS) [43]backbone optical communication networks. When tested on small networks against known results obtainable using unscalable methods like variants of integer/linear programming, it was shown to provide the optimal routing results.
The results provide the maximal number of communication pairs that could be accommodated given the network size, topology and number of wavelengths used; the minimal number of wavelengths required for a given network, topology and communication pairs; and the resulting utilization of edges. They identify the impact of using the suggested algorithm on the average path length, the utilization of wavelengths per edge and how it can be controlled using concave and convex objective functions.
In addition, our algorithms can be used for routing transmissions across networks in single instances, study the limitations of heterogeneous networks of different degree distributions, with variable edge signal-to-noise ratios and wavelength availability. Moreover, our algorithms could be employed in the design of new infrastructure, especially through the use of concave cost functions that consolidate routes, by determining the least important routes that could be removed with little effect on the network throughput or resilience. These are of both academic and practical values since the performance of optical communication networks is often directly related to their capacity limits, traffic congestion, rate of information flow and bandwidth flexibility. Moreover, we also studied a switching model, where wavelength can be converted (switched) at the vertices (transceivers) to settle an open and much debated question on the merit of wavelength-converters [44][45][46] for increasing throughput, resource utilization and resilience in optical communication networks.
While we mainly focus here on the optical communication network application and test the efficacy of the method on networks and number of wavelengths that are relevant to this application domain, one should point out that these problems are highly relevant to other domains. For instance, multi-wavelength NDP/EDP are relevant to both 5G and the future 6G wireless communication systems and wireless ad-hoc communication networks in the relay setting, where each node can act as a relay. Our algorithms can reduce path overlaps, which represent transmissions in similar wavelength, resulting in signal interference and low transmission quality, or to consolidate paths since longer paths result in signal degradation and the need for higher transmission power [26,27,47]. Another application is the design of very large system multilayer integrated circuits (VLSI), where non-overlapping wired paths to connect different components are sought to avoid cross-path interference.
In all cases, higher throughput, robustness and lower latency can be achieved for the same resource by obtaining a good approximation to the globally optimal solution. Practical algorithms for various applications often depend on the specific network topologies considered [48] and are aimed at maximizing the number of paths routed [49].
The reminder of the paper is organized as follows: in Sec. II we introduce the model used followed by the message-passing based algorithmic solutions for NDP, EDP and wavelength switching scenarios in Sec. III. Results obtained from numerical studies on a range of synthetic and real networks and a variety of objective functions are presented in Sec. IV followed by a discussion on their computational complexity. Possible extensions of the framework to accommodate real-world scenarios such as edges with different signal-to-noise ratios or wavelength availability are presented in Sec. V. Finally, we discuss the efficacy of the methods developed and point to future research directions in Sec. VI.

II. MODEL
We consider a dense wavelength-division multiplexing (DWDM) optical network G(V, E), with V ≡ {i | i ∈ G} and |V | = N , the set of nodes representing transceivers and E ≡ {(i, j) | (i, j) ∈ G} the set of edges, such that the indices (i, j) represent the optical fiber between node i and j. For a network which uses Q wavelength channels to deliver M transmissions, we introduce a variable s i,j on the link from node i to j such that s a i,j = s or −s if transmission s passes from node i to j or from j to i respectively, through link (i, j) using wavelength a; the transmission s corresponds to one of the origindestination pair {0, 1, . . . , M }. A similar variable s a i is defined for node i in the case of NDP, as explained in Sec. III. Since each transmission has to occupy an individual wavelength channel on a link, more wavelengths are generally required for more transmissions, i.e. a larger value of Q is required for a larger M . For specific network instances with Q wavelength channels, there exists a maximum number of transmissions denoted as M max which can be transmitted; alternatively, one can define the minimum number of wavelength channels, i.e. Q min , which accommodate all M transmissions on a specific instance. The relationship between M max and Q, or between M and Q min , would be highly relevant for characterizing the maximum capacity of optical networks.
This framework can accommodate a variety of objective functions; here, we consider the sum of the cost (or utility) on each link to be the objective function for optimization, given by where δ y x is the Kronecker delta such that δ y x = 1 if x = y and δ y x = 0 otherwise; the function F i,j denotes the cost on link (i, j) as a function of the argument in parenthesis; while it can take an arbitrary form we will mostly focus on simple polynomial functions. In the context of statistical physics, we introduce the inverse temperature β, and the partition function Z of the system is given by where Ω( s) is an indicator function such that Ω( s) = 1 if s satisfies all the constraints of the problem or otherwise Ω( s) = 0. The double vector notation comes to emphasize dependence on both topology and wavelength.
We now summarize the constraints of the optimization problem. Firstly, the route for each transmission must be contiguous. A loopless path is a sequence of nonrepeating nodes from origin to destination, for example, path s is constructed as O s → · · · → j → i → k → · · · → D s , where O s and D s are the origin and destination pair of transmission s. Secondly, for an intermediate node i along the path there exist only two used edges (j, i) and (i, k) for that wavelength and transmission. This constraint could be expressed as follows: if s a i,j = s = 0, then where ∂i ≡ {j | (i, j) ∈ E} is the set of the nearest neighbors of node i, and ∂i \ j ≡ ∂i − {j} is the subset of ∂i except node j. Regarding the wavelength constraints along the path, we will consider three scenarios in the subsequent analyses: 1. Node-disjoint paths (NDP) -where only a single transmission is allowed to utilize a specific wavelength channel on a node [8], but there can be multiple transmissions using different channels through the same node as shown in Fig. 1(a); this may correspond to transceivers which can only process a single transmission for each individual wavelength. The expression of the node-disjoint constraint is s a i,j = s = 0 ⇒ ∀k ∈ ∂i \ j : s a i,k ∈ {0, −s}. 2. Wavelength-switching (WS) with NDP -where all incoming or outgoing transmissions to transceivers (nodes) use different wavelength channels, but transmissions are allowed to switch between wavelength channels at the transceivers as shown in Fig. 1(b), leading to a larger routing flexibility.
3. Edge-disjoint paths (EDP) -where multiple transmissions using the same wavelength channel, are allowed to be routed through any given node [9]but cannot share an edge, as shown in Fig. 1(c). This is the typical scenario in optical communication networks.

III. MESSAGE-PASSING ALGORITHMS
To derive optimization algorithms to allocate simultaneously the optimal route and wavelength for a large number of transmissions, we solve the problem on multilayer graphs where each layer represents a different wavelength, and messages are passed within each layer for assignment of routes and between layers for allocation of wavelength. By applying the cavity approach from the study of spin glass systems [36,50], we can derive distributed message-passing algorithms for optimizing transmission routes in optical networks for NDP, WS and EDP scenarios.

A. Node-disjoint routing
First, we consider the NDP routing scenario [8] shown in Fig. 1(a). The network is represented by a factor graph with two types of variables, namely s a i defined on nodes and s a i,j defined on links. According to the Bethe-Peierls approximation [51,52], we assume only large loops exist in the network such that all neighboring nodes and edges of a node i are nearly independent in the absence of i.
To derive the multi-wavelength routing algorithm, we first define p i→j to be the message from node i to edge (i, j) and q i→j to be the message from edge (i, j) to node j. Both messages are conditional probabilities: p i→j (s) is the probability of edge (i, j) to be in state s due to the state of node i and q i→j (s) the probability of edge (i, j) to be in state s without the interaction from node i. We then write a closed set of recursion relations to express the message p i→j in terms of q k→i , as well as another set of relations representing q i→j in terms of p i→j and p j→i : where we have further defined the auxiliary quantities q b i,j given by For brevity we omit the normalization term and use the notation "∼" instead of the equality symbol. For a brief explanation of how these equations are construed: the first equation of (4), looks at the probability for no transmission on the edge i → j, as a summation of the probability of no transmission entering i (first term) and the probability of transmissions arriving at i from node m but leaving through some other node n (second term); the second equation looks at the message s arriving at node i but not leaving through any other edge but i → j; the third and fourth equations describe the probability of edge i → j to be in a given state 0 or s, given the related cost on the edge, in conjunction with all other wavelengths. In Eq. (4), node messages using wavelength a depend only on the messages from their neighbors using the same wavelength and the network is effectively mapped to a system with Q separate layers, each of which employs a different wavelength channel, as shown in Fig. 1(a); the interdependence between wavelengths is considered at the origin and destination of individual transmissions ( Fig. 2), as discussed below.
With the messages p i→j and q i→j having converged to stable values, one can express the marginal probability of node i being in state s using wavelength a as The marginal probability of edge (i, j) being in state s with wavelength a is given by where the factors after the dot symbols correspond to contributions from all neighboring links of (i, j) to the partition function, which takes into account all possible variable configurations in the trees terminated at the neighboring edges.
To simplify the algorithms, we introduce the variables φ ≡ − 1 β log q and ψ ≡ − 1 β log p, and take the zerotemperature limit for optimization β → ∞. The recursion relations of messages in Eq. (4) thus become equivalent to the min-sum belief propagation relations, given For each individual transmission, we then introduce auxiliary nodes labeled as µ = ±|µ|, which connect to the origin and destination nodes in each of the Q wavelength network layers, respectively (see in Fig. 2). These auxiliary origin-destination pairs of each transmission communicate with all wavelength replica networks to determine the path and wavelength channel allocated to each transmission, using the message passing algorithm, given by such that the transmission routes are determined independently on each wavelength channel. After introducing the auxiliary nodes, we treat the messages to and from them as to other network nodes. For example, the neighboring nodes set of node i in Fig. 2 is ∂i = {j, k, l, µ}, and the calculation of messages i → j and i → µ both follow Eq. (8). The marginal probability of edge (i, j) being in state s using wavelength a is then given by and the state s a i,j of wavelength a of edge (i, j) is deter-mined by which ultimately leads to the optimized configuration of routes for all M transmissions through the optical network with Q wavelength channels.

Linear cost
While the objective function can take many different forms we use a simple power as the cost function F i,j (x) on edges x γ [4,53], as it provides a good example of both concave and convex costs. When γ = 1, the cost function is linear and equivalent to the total length of all transmissions L ≡ (i,j);a 1 − δ 0 s a i,j , and optimizing it is equivalent to finding the shortest average path.
In this case (γ = 1), the message-passing equations Eq. (8) can be simplified to where and the variable ψ in Eq. (8) can be omitted. The marginal quantities are then given by and the messages to the auxiliary nodes at the origins and destinations are the same as Eq. (10).

Switching wavelength channels at nodes
A generalization of the NDP scenario is to allow for transmissions to change wavelengths at nodes (WS) as shown in Fig. 1(b), which is equivalent to a singlewavelength network routing with both node and edge capacity being Q, and utilize other multiplexing techniques than wavelength division (e.g. code or time division multiplexing). A variable τ is introduced to enforce this capacity constraint on nodes. In this case, the message-passing equations are given by where τ µ i denotes the summation overψ i from a total of at most Q transmissions, excluding µ, that pass through node i, given by Specifically, τ µ i (1) could be understood as indicating that there are free wavelength channels on i that µ could take, and τ µ i (0) indicates that transmission µ does not pass through node i. The quantityψ µ i is an auxiliary variable related to the probability that node i would be a part of the path of transmission µ, which leads to The messages from the origin and the destination of transmission µ to the neighboring edges are slightly different from the previous case of NDP without switching. In Eq. (16) Q i = Q, andψ i (1) = 0 andψ(0) = ∞ for the transmissions using node i as their origin or destination. Alternatively, one can set Q i = Q − M i , where M i is the number of transmissions using node i as their origin or destination, and then exclude these transmissions when one calculates τ i . Here, we first introduce a variable µ i = ±1 to denote node i as the origin or destination of transmission µ respectively, given by The variable φ µ i,j determines the state of edge (i, j), i.e. whether it is a part of the path of transmission µ or not, and its expression is given by such that the optimized state is σ µ i,j = arg min σ φ µ i,j (σ).

B. Edge-disjoint routing
Edge-disjoint path routing (EDP) is similar to NDP routing but is less restrictive, since nodes can accommodate any number of paths with the same wavelength but edges do not [9] (Fig. 1(c)). In other words, for node i, if there exists an edge (i, j) with s a i,j = s 0 = 0, then there exists one and only one edge (i, k) with s a i,k = −s 0 continuity of transmission path, while other neighboring edges using the same wavelength channel a could be either in state 0 or take up other transmissions such that s a i,m1 = −s a i,n1 = s 1 = 0, s a i,m2 = −s a i,n2 = s 2 = 0 etc. In comparison, for NDP at most two variables s a i,j for neighboring nodes j ∈ ∂i can assume a non-zero value. In other words, single-wavelength EDP is equivalent to a generalized version of WS, where the capacity of each node is non-uniform and determined by its degree and the number of it being chosen as origins or destinations of transmissions. Consequently, the message passing equations for EDP scenarios are more complicated and their computational complexity is higher. The corresponding message passing equations are given by where the message from node i to edge (i, j) (first two equations) are different from those of NDP in Eq. (8), while Eq. (9)- (12) in the NDP scenario apply also here. It is numerically difficult to compute the message passing relations in Eq. (20), since all feasible configurations satisfying the constraints have to be considered by the first two equations, which results in a computational complexity of where K = |∂i| − 1 is the number of neighboring edges of node i, M is the total number of transmissions, x is the floor function that is equal to the greatest integer less than or equal to x, n P k = n! (n−k)! is the number of ordered permutations of k out of n elements and n C k = n! k!(n−k)! the number of unordered combinations. To simplify the computation we map the task on to the maximum weight matching problem [9,54]. As in [9], we consider all the possible pairs of considered edges {(i, k) | k ∈ ∂i \ j} and obtain the weight of each pair Then, we construct a weighted graph with |∂i| − 1 nodes where the weight of each edge is given by Eq. (22). After that, ψ a i→j in Eq. (20) could be obtained by the maximum weight matching algorithm [9,55].
In Fig. 3, we show the transformation employed for calculating the messages i → j. Links between all pairs, including auxiliary origin-destination and graph edges, are shown on the left. The right figure shows all valid links between nodes, while distinguishing between auxiliary and graph nodes: any of the ordinary graph nodes could be paired but the auxiliary nodes could only be paired with ordinary graph nodes. The maximum matching algorithm finds edge sets with the maximum sum of weights, where matched edges have no common nodes; the inverse value of the weights sum is approximately the value of minimum matched configurations in Eq. (20).

Linear cost
Similar to Sec. III A 1, if the cost function on edges is linear F i,j (x) = x, then the message-passing equations could be simplified as follows, The marginal states of nodes and edges, and messages to the auxiliary nodes at transmission origins and destinations, are the same as in the NDP case of Eq. (10) and (14).

C. Algorithmic framework
The same algorithmic procedure governs both NDP and EDP routing algorithms based on the message passing equations. The min-sum algorithm, i.e. zerotemperature optimization algorithm, follows the process outlined below: 1. initialize the messages {ψ a i→j , φ a i→j , φ µ→a }; 2. update the messages by Eq. (8)-(10) in the nodedisjoint scenarios, and Eq. (20) and (10) in the edge-disjoint scenarios until convergence or when a maximum number of iteration steps is reached; 3. calculate the marginal state of edges by Eq. (11) and the state of each wavelength channel by Eq. (12).
After obtaining the marginal states of edges and wavelengths, the allocation of transmissions to specific wavelength and path would follow. However, there are some cases when the algorithms need many iterations to provide a valid configuration, for example if the cost function F i,j (x) in Eq. (1) is concave, e.g. F (x) = √ x in these cases, a decimation procedure can be introduced to speed up convergence. As an example, here we show the NDP algorithm with decimation (fixing states at intermediate steps) incorporated as part of the process: 1. initialize the messages {ψ a i→j , φ a i→j , φ µ→a } and the state of all wavelength channels per available (undecimated) edge; 2. update the messages of the available wavelength channels and edges by Eq. (8)-(10) for a specific number of iteration steps; 3. compute the marginal states of the available wavelength channels and edges by Eq. (11); calculate the quantityφ a i,j = min s =0 φ a i,j (s) − φ a i,j (0), and fix the state of the wavelength channel a of edge (i, j) with the largest value ofφ a i,j to be s a i,j = 0.
4. iterate step 2 and 3 and determine the state of the available wavelength channels and edges by Eq. (12) until a valid solution is obtained or a maximum number of iteration steps is reached.
The algorithms with or without the decimation procedures have a similar performance in terms of the optimized cost, but for harder problem (e.g. F (x) = √ x) the decimation procedures can reduce the number of iteration steps needed for generating a valid routing solution.

A. Linear cost
We first examine the NDP, WS and EDP scenarios with linear cost F (x) = x to explore the behavior and performance of communication networks, such as capacity and average length, with impact on latency and number of wavelength channels required.

Dependence on the number of wavelength channels
We performed numerical experiments using the three algorithms on three different types of generated random networks including random regular graphs, Erdős-Rényi and scale-free networks with 100 nodes and an average node degree of 3, as well as on two real optical communication networks -CONUS60, which has 60 nodes and 79 edges [43] and BT-Core, 22 nodes and 35 edges [56] as shown in Fig. 4. Specifically, for the scale-free networks studied, the degree distribution is p(d) ∼ d −2.63 .  We define the capacity to be the maximum number of transmissions, which we denote as M max , that can be transmitted by an optical communication network with Q wavelength channels. We remark that capacity depends on network topologies and the set of origin-destination pairs, so in Fig. 5 we report the dependence of the average value of M max with standard deviation as error-bars on Q, showing the average behavior of capacities for different Q values.   5 shows the simulation results for the five types of networks -in the three different routing scenarios. For the three types of random networks and BT-Core, the average capacity per wavelength channel M max /Q keeps increasing at the beginning and become saturated as the number of wavelength channels Q increases.As for the networks CONUS and BT-Core in Fig. 5, we see that the average capacities also increase with Q but the increases are not as fast as that in the three generated networks.
The average path length of the corresponding networks are shown in Fig. 6 and 7. In Fig. 6, we see that as more wavelengths are available the shorter the average path length L/M decreases; for instance, in NDP scenarios shown in Fig. 6(a)-6(c), the average path length with Q = 1 is higher than that with Q = 4, which is higher than that of Q = 9 as the average load per wavelength channel (M/Q) increases. Moreover, with the same values of Q and M , it is clear that the wavelengthswitching model would provide us with shorter paths and a higher capacity on Erdős-Rényi and scale-free networks as shown in Fig. 6.
In comparison with NDP, the advantages of EDP with multiple transmissions using the same wavelength availability is clearer on random regular graphs as shown in Fig. 7(a), and less so on other graph types, presumably due to the graph heterogeneity and finite size effects. We observe that EDP routing increases the network capacity, i.e. average transmission load per wavelength channel, and provides shorter valid paths.

Dependence on the number of transmissions
In practice, one may wish to minimize the total number of active wavelength channels for given demand. We define the smallest number of wavelength channels which accommodate a set of transmission to be Q min ; to find Q min from for a specific instance, we gradually increase Q  We compare the values of Q min /M obtained by our algorithms and by the MGA algorithm on random regular networks, the CONUS and BT-Core networks, in both NDP and EDP scenarios. The numerical results on random regular networks presented in Fig. 8(a) and 8(c) show that our algorithms offer significant advantages over the MGA algorithm in reducing the number of wavelength channels required for specific random transmission pairs in. The improved performance is also observed on the CONUS and BT-Core networks in Fig. 9 and 10 respectively. The underlying reason for the improvement is the low success rate of random greedy assignments. The numerical results demonstrate that our MWR algorithms lead to a better use of resource. The additional flexibility provided by the wavelength-switching on transceivers leads to an even smaller value of Q min as M increases in Fig. 8(a) and 9(a), whereas the experiments on the BT-Core network do not show a significant improvement of WS over the ordinary NDP in Fig. 10(a), which may depend on the topology of the graph or our range of values tested.
Comparing the average path length found for both NDP and EDP routing on random regular graphs, we see in Fig. 8(b) that the average path lengths obtained by the MGA and MWR for NDP routing are similar but that shorter paths are found by MWR in the EDP scenario ( Fig. 8(d)). A significant reduction in average route length is not expected since typical route lengths on random graphss is O(log N ) with little variability. However, a significant reduction in the number of wavelength channels used Q min , is shown in Fig. 8(d) and a similar trend is shown for EDP scenarios in Fig. 9(a,c) and 10(a,c) for the two real networks.

B. Computational complexity
Next, we discuss the computational complexity of our multi-wavelength routing algorithms with a linear cost function F i,j (x) = x. In this case, if we use the normalization φ a i→j (s) := φ a i→j (s) − φ a i→j (0) in the min-sum equations (13) (15) (20), the computational complexities would decrease.
Under NDP routing, the complexity in computing Eq. (13) is approximately O( k 2 M ), where k is the average node degree, and there are k N Q such messages. The complexity in computing the messages from origin/destination nodes in (10)  , which linearly scales with Q as shown in Fig. 11(a).
For EDP routing, the complexity of matching a pair of incoming and outgoing transmissions with the same wavelength in Eq.(23) is O( k 3 log k ) [9]. Therefore, the total complexity of one round of update for all messages is O(2M k 4 log k × k N Q + M Q 2 ), and if the graph is sparse the complexity becomes O(M N Q + M Q 2 ). For relatively small Q, the complexity scales with Q as in NDP routing as shown in Fig. 11(b).
However, when M/N N , the effective degree becomes quite different from the average degree k ⇒ k + 2M N (see in Fig. 2, the degree of node i increases from 3 to 4 after introducing the auxiliary node µ). When M is very large, for instance M ∼ O(N 2 ) as is used in many optical communication network applications, the effect of effective degree cannot be omitted. To make the algorithm scale better we divide the messages into three types: 1. messages from auxiliary nodes to ordinary nodes, e.g. µ → i in Fig. 2 -there are 2M Q such messages; 2. messages from ordinary nodes to ordinary nodes, e.g. i → j in Fig. 2 -there are k N Q messages of this type; 3. messages from ordinary nodes to auxiliary nodes, e.g. i → µ in Fig. 2 -there are also 2M Q such entities.
The messages of type 1 obeys Eq. (10) and the resulting total complexity is O(M Q 2 ). Noticing that the messages from auxiliary nodes are sparse, only φ µ→a (0) and φ µ→a (µ) are non-trivial in Eq. (10), and only φ b→µ (−µ) is needed to generate new messages, which simplifies the calculations of Eq. (13). Under NDP routing, by separating messages of type 1 and 2, the computational complexity of generating one message of type 2 or 3 by Eq. (13) is O( k 2 M ). Then, the total complexity of the three types of messages is  (N 5 ). In Fig. 11(c), we present the one-iteration computing time on random regular graphs with degree k = 3 and N nodes, the one-iteration computing time scales approximately as N 5.14 , which fits our analysis of Under EDP routing, the computational complexity of type 2 and 3 messages depends on the complexity of the maximum weighted matching algorithm (approximately O(( k + 2M N ) α )) with α being a coefficient to be determined. Using the sparse properties of type 1 and 3 messages, the total complexity of the three types is O(M Q(N + Q + (M/N ) α )). Recalling the mapping to matching problems -the matching network has approximately k + 2M N nodes and k ( k −1 2 + 2M N ) edges, and the average node degree is approximately k when 2M N k , therefore it is a sparse network. By using the algorithm in [9], the complexity should be approxi- ). However, due to resulting topology and the special network structure our tests result in complexity scaling close to O( k + 2M N ) when tested. In other words, α = 1 and therefore the resulting total complexity is O(M Q(N + Q + M/N )). Experiments on random regular graphs presented in Fig. 11

C. Non-linear cost
In this section, we show the simulation results where a more general form of the cost function where the argument is the wavelength occupancy q i,j on edge (i, j), given by When γ > 1, the utility increases faster on heavily loaded edges, biasing solutions towards paths with uniform load, balancing the loads on edges. On the other hand, when 0 < γ < 1, the cost on edges increases slower with the load, leading to configurations which consolidate transmissions on used edges, leaving more edges and wavelength channels idle. Such configurations would be relevant for identifying less important transceiver nodes which can be switched off in hours of low usage. In Figure 12 we show the distribution p(x) of wavelength occupancy on edges of random regular graphs and the two real network CONUS and BT-Core for different values of γ. For all three networks and both NDP and EDP scenarios, we observe that when γ = 0.5 more edges were unused, i.e. a higher value at p(x = 0) in Fig. 12, and as γ increases the distribution become more evenly distributed and peaked at some values of x, which corresponds to the balancing of edge loads. We show in the insets the relative difference between p(x) obtained by γ = 1 and those obtained by γ = 0.5 or 2, i.e. δ(x; γ) ≡ p(x;γ)−p(x;γ=1) p(x;γ=1) , where the evidence for load consolidation or balancing, for γ = 0.5 or 2 respectively, become obvious.

D. Comparison with linear programming
Linear programming is a commonly used method for optimized throughput of optical networks [57]. Here, we compare the results of our proposed algorithm with those To solve the routing problem using linear programming, we find n shortest paths for each of the M transmissions, denoted by the variables { s µ k } with superscript µ = 1, . . . , M representing transmissions and subscript k = 1, . . . , n representing the k th candidate path for transmission µ. By carrying either a subscript for nodes or edges, the variables s µ k can represent a configuration of node states or edge states, and the NDP or EDP constraints can both be expressed in terms of s µ k . In this case, we introduce a variable σ a,µ k = 1 to denote that transmission µ chooses its k th path with wavelength a, σ a,µ k = 0 otherwise. The disjoint constraint can be ex-pressed as follows, In Eq. (25), if the path s µ k represents a configuration of node states, then each element of v a represents the load of a th wavelength channel on that node, and node-disjoint constraint restricts the value of the load to be no more than 1; if the path s µ k represents a configuration of edges, then the edge-disjoint constraints can be defined in a similar manner.
For each transmission, we choose one wavelength to accommodate one candidate path, which is given by The objective function for the problem with linear cost is given by The problem can be solved by linear programming with the objective function of Eq. (27) subject to the constraints in Eq. (25)- (26) where all the expressions are linear.
For the WS scenario, the variables σ µ k are introduced instead of the previous σ a,µ k , and the capacity constraint for each node is given by where s µ k is a configuration of node states. The constraint for all the transmissions is given by The objective function to be minimized is j v j . We conducted numerical experiments on four real networks [56,57] and compared the results obtained by our algorithms and linear programming in the NDP, wavelength-switching and EDP scenarios. We show the smallest numbers of wavelength channels, i.e. Q min , needed to complete all possible transmissions, and the corresponding total path length L in Tab. I. As we can see, our message-passing algorithms and linear programming yield almost identical performance in finding optimized path solutions. For relatively small network, linear programming is efficient, but it quickly becomes impractical when the size of networks increases, whereas we show in Sec. IV B that our message-passing algorithms have more practical scaling properties. Table I. The smallest number of wavelength channels Qmin required to transmit all M = |V |(|V | − 1)/2(= N (N − 1)/2) transmissions on four small real optical communication networks including NSF-Net, Google-B4, DTAG/T-systems and BT-Core, obtained by our multi-wavelength routing algorithm (MP) compared to that obtained by linear programming (LP), in edge-disjoint (ED), node-disjoint (ND) and node-disjoint wavelength-switching (WS) scenarios. The corresponding total path lengths are also shown. In real optical networks, the number of wavelength channels in different optical fibers may vary; their lengths, signal-to-noise ratios, or type of fibers used may influence the quality of communication. To make our model more general and realistic, we consider the case of optical networks with weight w i,j and number of wavelengths Q i,j defined for any individual link (i, j) (or Q i defined for node i in node-disjoint wavelength-switching scenarios). With simple modifications, our proposed algorithm can be generalized to accommodate w i,j , Q i,j or Q i .
In cases with heterogeneous Q i,j on edges, we denote the largest number of wavelength channels among all edges to be Q * , i.e. Q * = max (i,j) Q i,j , then for an edge (i, j), one can introduce Q * − Q i,j additional wavelength channels with state 0, such that all edges on the network would have virtually Q * wavelengths.
As for heterogeneous weights on edges with linear cost, the objective function becomes L = (i,j) w i,j A. Heterogeneous wavelength availability on edges In the message-passing equations (13) describing NDP scenario, the first equation φ a i→j (0) considers the condition for wavelength a of edge (i, j) to be in state 0 in the absence of node j, and the contribution to the objective function L is 0; this equation is the same in cases with heterogeneous edge weights. Nevertheless, weights w i,j should be introduced in the second equation as The marginal messages on edges, Eq. (14), should be modified as φ a i,j (s) = φ a i→j (s) + φ a j→i (−s) + w i,j (δ 0 s − 1).
The same algorithmic procedure described in Sec. III C can be applied in the present cases with heterogeneous edge weights.
For an even more general scenarios, our model can be modified to accommodate heterogeneous weights for different wavelength channels on the same edge, i.e. w a i,j = w b i,j for a = b, asymmetric directed weights such as w i,j = w j,i on directed graphs. In these cases, one only needs to replace w i,j in the above modified equations by w a i,j or directed weights.

B. Heterogeneous wavelength availability on nodes with wavelength-switching
For the NDP scenarios with heterogeneous wavelength availability on nodes, we do not have to introduce additional wavelengths and keep unavailable channels in state 0 if Q i are non-uniform; this is because the equation of node capacity constraint Eq. (16) has already considered the case of different Q i values. Only the second equation of Eq. (15) needs to be modified, which reads φ µ i→j (±1) ∼ w i,j + τ µ i (1)+ min

VI. CONCLUSION
Multi-wavelength NDP/EDP routing lies at the heart of the efficient running and design of optical communi-cation networks, that act as the backbone of the Internet. One of the key questions in running optical communication networks more efficiently is in the ability to carry out these routing tasks effectively for large systems. This serves for both day-to-day running of the network as well as for the design of new networks and the modification of existing infrastructure. While principled single wavelength NDP and EDP routing algorithms based on message passing have been developed already, they could not be employed in real optical networks due to the difficulty in extending the algorithms from the single wavelength to the multi-wavelength case. This essential and pivotal aspect of routing in optical communication networks makes the existing single-wavelength methods intractable and requires computational cost that grows exponentially with the system size.
To accommodate a large number of wavelengths and transmissions in large systems, we have developed algorithmic solutions that include multi-layer graphs, where each layer represents a different wavelength, and messages are passed within layer (routing assignment) and between layers (wavelength allocation). The scalable algorithm we have devised shows very good performance in manageable time scales.
We expect the algorithm to be implemented in realistic scenarios, where specific aspects of real network routing, such as heterogeneous wavelength availability and signal-to-noise ratios will have to be added in the man-ner outlined in Sec. V. We also expect the algorithm to be utilized for network design and see several possible extensions for both localized and global message passing-based implementation. Utilization of the algorithms developed here in ad-hoc network communication, multilayer VLSI design and multilayer networks will require further study of the specific requirements for the different applications.