Search Trajectory Networks Applied to the Cyclic Bandwidth Sum Problem

Search trajectory networks (STNs) were proposed as a tool to analyze the behavior of metaheuristics in relation to their exploration ability and the search space regions they traverse. The technique derives from the study of fitness landscapes using local optima networks (LONs). STNs are related to LONs in that both are built as graphs, modelling the transitions among solutions or group of solutions in the search space. The key difference is that STN nodes can represent solutions or groups of solutions that are not necessarily locally optimal. This work presents an STN-based study for a particular combinatorial optimization problem, the cyclic bandwidth sum minimization. STNs were employed to analyze the two leading algorithms for this problem: a memetic algorithm and a hyperheuristic memetic algorithm. We also propose a novel grouping method for STNs that can be generally applied to both continuous and combinatorial spaces.


I. INTRODUCTION
O BSERVING the inner dynamics of metaheuristics is crucial for a better understanding of how these algorithms differentiate from each other on the way they explore the search space, and how those differences relate to their performance under specific scenarios. This is increasing in relevance, as in recent years numerous novel metaheuristics have been proposed. The metaphorical formulation of these algorithms often takes inspiration from processes found in natural systems, physics [1]- [3], social behaviors of species like ants [4], bees [5], fireflies [6], chaos and game theory [7]- [9], etc. Furthermore, there is a tendency for metaheuristic hybridization [10]- [12] and hyperheuristics [13], [14]. Very promising results have been obtained by such proposals, causing a growing interest in extending their applications and also increasing the need for analysis tools to effectively characterize their behavior.
On this matter, search trajectory networks (STNs) [15], [16] are a relatively novel analysis tool for gaining a better perspective on the inner search dynamics of metaheuristics in relation to how they explore the search space of a particular problem instance. STNs allow to identify how the success of a metaheuristic relates to the search process being conducted through specific search regions. Therefore, STNs constitute a metaphor-free attempt at profiling metaheuristics.
Both local optima networks (LONs) [15], [17] and STNs represent, in the form of directed graphs, key structural features of the search space and how the algorithms navigate them. Both models are built as networks and analyzed using network metrics to identify the paths conducting towards the best search space regions, which areas are hard to escape from, and how they are related. The key difference is that a search trajectory network is created from trajectories between solutions that are not necessarily locally optimal. In that way, STNs overcome some limitations of LONs, opening their applicability to contrast the behavior of a wider variety of metaheuristics, such as population metaheuristics that do not necessarily perform local search.
STNs were initially proposed to characterize differential evolution and particle swarm optimization [18] for several classical continuous optimization benchmark functions [19]- [21]. In this scenario the continuous search space was partitioned into regular size hypercubes. STN analysis was later extended to cover not only population-based algorithms but also stochastic local search methods, and both continuous and combinatorial optimization problems [16]. However, representative standard metaheuristics, as opposed to stateof-the-art methods, were considered. As the combinatorial case study, a classic facility location problem [22] (p-median) encoded with binary strings was used. In the continuous domain, search can be partitioned into discrete hypercubes of a predefined length, thus defining the STN locations. This cannot be done for discrete spaces. For discrete spaces, locations can be associated with individual solutions. However, it is still desirable to coarsen the search space in order to appreciate the trajectories at different levels of granularity. Therefore, the authors proposed a partitioning scheme based on the sampled solutions [16]. The idea was to group as a single location all solutions sharing the same value in some variables while ignoring (removing) other variables. This method worked well for a binary representation, but it is not suitable for a permutation encoding.
In this work we report a STN-based analysis of a combinatorial optimization problem, the cyclic bandwidth sum problem (CBSP) [23]. Our goal is to gain a new perspective on how the two state-of-the-art metaheuristics proposed for the CBSP behave, why are they successful, why some type of instances seem more challenging, as well as obtaining new insights of the underlying fitness landscapes. These algorithms are interesting subjects for STN analysis because their design features differ from the classical memetic algorithm they were originally based on. The first algorithm is a memetic algorithm (MA) that employs a nonintensive local search approach [24], the second one is a hyperheuristic approach [25] using a dynamic multi-armed bandit (DMAB) [26] to automate the selection of the genetic operators and the fitness function [27] within a MA.
Another contribution of this work is the proposal of a search space partitioning method based on reducing the solutions dimensions while preserving their distances using multidimensional scaling [28]- [30]. This approach to search space partitioning is more general than those previously used for STN analysis [15], [16], as it can be applied to any form of solution encoding (discrete and continuous) as soon as a suitable distance metric between solutions can be devised.
The remaining sections are organized as follows. Section II provides a brief introduction to graph embedding problems and in particular to the CBSP. Section III describes the methodology for the creation and analysis of STNs. Section IV introduces the two studied algorithms. The experimental results are presented and discussed in Section V. Finally, Section VI presents our conclusions.

II. THE CYCLIC BANDWIDTH SUM PROBLEM
Graph embedding problems (GEPs) [31] are a family of combinatorial optimization problems focused on the rearrangement of graphs to fit a certain new layout. These type of problems are commonly defined in terms of a guest graph and a host graph. The guest graph describes the original relationships among whichever objects the graph represents, for example, a virtual computer network or a set of related facilities. Meanwhile, the host graph describes the new layout to rearrange the guest graph, for example, a network infrastructure architecture, or a set of connected physical locations for facilities. Therefore, graph embeddings are equivalent to mappings that assign guest nodes to host nodes, and guest edges to paths in the host. The embeddings are often represented as labelings [32], assigning to each guest node a host node represented by using its label.
Commonly, a GEP's objective is to find an embedding that optimizes a certain measure, defined with respect to the topological structure of the host graph. The generalized classification of GEPs [32] includes three groups characterized by the nature of what are they trying to optimize: a) the distance among guest nodes when embedded in host nodes, b) the sum of distances among guest nodes when embedded in host nodes, and c) the cutting of host edges embedded in host paths.
One of the best-known and widely studied GEP is the bandwidth problem [41], consisting in rearranging a graph into a linear layout, i.e, embedding a guest graph of order n into the path graph P n in such a way that the linear distance among adjacent guest nodes, i.e., the bandwidth, is minimized. The bandwidth problem belongs to the first group mentioned above. The bandwidth sum problem [42]- [44] is then the GEP consisting on minimizing the sum of the linear distances. Now, consider how a cut on a particular edge of the host graph will disrupt the embedded guest edges that cross through it. The cutwidth problem [32], [45], [46] consists in minimizing the number of disrupted guest edges affected by a cut on any host edge.
Particular GEPs are also defined by the topology of the host graph. The examples mentioned above involve only linear topologies in the form of path graphs as hosts. When the guest graph topology is instead circular, denoted by the cycle graph, the arising problems are the cyclic bandwidth [47], the cyclic bandwidth sum [48], and the cyclic cutwidth [49], respectively. This paper deals with the cyclic bandwidth sum problem (CBSP) [48], where the objective is to find the optimal way to embed a graph into a circular layout while minimizing the sum of cyclic distances between adjacent guest nodes.
Formally, the CBSP is defined as follows. Let G = (V, E) be a finite undirected graph (the guest) of order n = |V | and C n a cycle graph (the host) with vertex set |V H | = n and edge set E H . Given an injection ϕ : V → V H , representing an embedding of G into C n , the cyclic bandwidth sum (the cost) for G with respect to ϕ is defined as: where |x| n = min{ |x|, n − |x| } (with 1 ≤ |x| ≤ n − 1) is called the cyclic distance, and the label associated to vertex u is denoted ϕ(u). The CBSP consists in finding the optimal embedding ϕ * , such that Cbs(G, ϕ * ) is minimum, i.e., ϕ * = arg min ϕ∈Φ {Cbs(G, ϕ)} with Φ denoting the set of all the potential embeddings. The value of the optimum can be calculated in an exact way for the following graph topologies: path, cycle, complete bipartite graph, wheel, and power of cycles [23]. Approximations of the optimal value were reported for Cartesian products of two graphs (paths, cycles, and complete graphs) [50]. Regarding heuristic solution proposals, the problem has been approached by general variable neighborhood search [51], a constructive greedy heuristic [52], as well as the MA [24], and DMAB+MA [25] algorithms analyzed in this work.

III. BUILDING SEARCH TRAJECTORY NETWORKS
Given a metaheuristic algorithm A, an STN is built as a weighted, directed graph ST N A (V A , E A ) where: • V A is the set of search space locations visited by the algorithm. • E A is the set of directed edges, such that two nodes a, b ∈ V A are adjacent if the algorithm A performed a transition between solutions in the respective locations. The weight of the edge represents the number of times the transition from a to b occurred during the search. A location is a search space region that contains at least one representative solution. Each representative solution was the current best solution during an iteration of the algorithm, corresponding to the best solution in the population for multi-trajectory metaheuristics or the incumbent solution for single-trajectory metaheuristics.
In order to compare the search dynamics of two or more metaheuristics, the STN models for two or more metaheuristics, under the same problem instance can be merged. This idea of merging the trajectories of two algorithms in a single network model was first proposed and implemented for local optima networks in [53], and later adapted to STNs [15], [16].
More formally, the merged STN for two metaheuristic algorithms A i and A j is the union of the STN graphs ST N Ai The weight of each edge in the merged STN is equal to the sum of the number of times algorithms A i and A j performed the transition represented by the edge.

A. SAMPLING, MODELING AND PARTITIONING
The first step for creating an STN is to record an aggregated sample of the transitions between representative solutions that occurred during independent executions of the studied metaheuristic for a problem instance. A parameter is employed to control the frequency of recording those transitions. This process is quite simple to implement, it has no effect on the metaheuristic search process, and it does not represent a significant overhead.
The STN modeling is the process of turning the sampled trajectories into a network of locations, under the previous definition of nodes, edges and weights. An important part of this process is to establishing a mapping between single solutions and locations. By definition, a location can be a subregion as small as containing only one single solution. However, the STN can result more informative if locations represent instead subregions of the search space. The definition of which solutions are mapped to the same location can vary. In general, it should be linked to a notion of closeness among solutions. For example, for continuous optimization problems, a partitioning of the search space into hypercubes of regular size has been induced by adjusting the precision of the solution encoding [15].
Since the quality of each location is represented by the fitness of the best solution that belongs to it, it is implied that, if the metaheuristic can reach a solution within the corresponding subregion, it is likely that it could as well reach the best solution from it. Grouping solutions into STN locations also provides a globalized perspective on the overall search dynamics. For example, by making it possible to identify if an algorithm is recurrently visiting some specific regions, and to detect regions that act like hubs, even if the algorithm does not revisit the exact same solutions.
The partitioning method employed in this work is based on the classical multidimensional scaling [28], [29]. The idea is to reduce a set of solutions from a representation in n dimensions to the Euclidean space, while minimizing the distortion of the pair-wise distances among them. The solutions we dealt with in both algorithms were encoded as permutations of n elements. Under the problem definition, they are circular permutations. Therefore, we employed the interchange distance [54] to evaluate the pair-wise distances between unique solutions in the sample. The interchange distance measures the number of cycles between two permutations, therefore it is suitable to reflect the swap-based edition distance for CBSP potential solutions. Then, each solution was mapped to a 3-tuple (three dimensional coordinates) by using classical multidimensional scaling. At that moment, the obtained coordinate values are rounded up to integers for grouping solutions sharing the same resulting integer 3-tuples. Multidimensional scaling has been previously employed in the creation of visual representations of the distribution of solutions to analyze the cartography of search spaces [30], [55].

B. ANALYZING STNs
The analysis of STNs has two complementary components: a set of network metrics and a visual representation. The metrics capture relevant features of the merged networks and their structure, while the visual representation provides a better idea of how the trajectories intersect and allows to observe relevant locations.
The set of network metrics evaluated in this work are the following: • nodes: Number of nodes • edges: Number of edges VOLUME 4, 2016 • end nodes: The number of unique nodes without any edge going out from them. Each end node corresponds to the end of a trajectory. • best nodes: The number of unique nodes with a fitness equal to the best found during the trajectory. • in-strength: Evaluated for the best nodes and end nodes, it is the ratio between the sum of their weighted incoming degree and the sum of weighed incoming degrees for all end nodes. • shared nodes: Number of nodes visited by more than one algorithm. • visited nodes (v. nodes): Number of nodes visited by algorithm A i . • shortest path length (s. path): Length of the shortest path to best node visited by algorithm A i . • number of shortest paths (n. paths): Number of shortest paths to the closest best node visited by algorithm A i . The STNs were plotted in R by employing the igraph package [56] and the force-directed layout algorithms Kamada-Kawai [57] and Fruchterman-Reingold [58]. Both layout algorithms try to create displays of the network that have: • a roughly even distribution of vertices • minimized crossings of edges • approximately uniform length edges • preservation of inherent symmetries, in such a way that similar subnetworks are depicted similarly as well.

IV. THE HYBRID ALGORITHMS
The algorithms considered in this study are a memetic algorithm (MA) [24], and DMAB+MA [25], a technique that employs a dynamic multi-armed bandit [26] as a hyperheuristic approach to adaptive operator selection (AOS) [59] within a MA. DMAB+MA can be considered also as an adaptive memetic algorithm. Both are hybrid populationbased algorithms, and to the best of our knowledge, they are the top state-of-the-art methods regarding solution quality for the CBSP. These algorithms are the result of prior research on the interplay of different configurations of genetic operators within a MA, extending to the selection, crossover, fitness function, mutation and survival strategy, and how the strength of individual configurations can be combined in an AOS approach.

A. MEMETIC ALGORITHM
The MA structure is described in Algorithm 1. The MA begins with a population P of µ permutation encoded randomly generated parents. Throughout the generations, pairs of parents are selected using binary tournament. Each pair of parents produces one offspring individual o with probability p c , by using cyclic crossover. Alternatively, with probability 1 − p c , individual o is instead a copy of the fittest parent in the pair. Individual o is mutated with probability p m . The mutation operator is a fitness oriented reduced 3-swap that works as follows. Three random distinct genes are chosen, then the impact over the fitness of the six different ways on which those genes can be swapped in pairs are calculated. The mutated individual o results from performing the swap that result on the best fitness improvement, or the one that deteriorate the fitness the least (if there is not an improving swap). After this, an inversion operator is applied with an independent probability p i , producing a further mutated individual o . Together, the reduced 3-swap and insertion operators can be seen as a two phase mutation in which each phase is independent. The first phase is oriented to seek good fitness mutations, while the second has a random nature. While only individual o joins the offspring population O, possible improvements on the best-so-far solution achieved by individuals o and o are still recorded. This MA uses the (µ, λ) type survival, directly replacing the parent population P by the offspring population O.
The local search phase occurs after the survival. Under this MA approach, the local search is non-intensive, in the sense that a) it is not applied to each individual in the population, but rather only to the fittest individual in it, and b) it is not carried on until a local optimum is reached, but instead only for a reduced number ls of iterations. These particularities are intended to prevent the population from becoming rapidly overtaken by locally optimal individuals, whose genes could proliferate excessively and cause premature convergence. Notice that the local search can resume its process where it was left if the best individual in the population remains the same for more than one generation (i.e., P best,t = P best,t+1 ), but also it can be restarted from similar fitness value individuals (i.e., f (P best,t ) = f (P best,t+1 ) and P best,t = P best,t+1 ). P best ← localsearch(P best , ls) 18: g ← fitter individual among current g and P best 19: until stop criterion is met 20: return g

B. DMAB+MA
The multi-armed bandit (MAB) [60]- [62], is an effective way to approach the problem of choosing the best among k alternatives. The MAB is traditionally modeled as k casino bandit machines, each with an underlying unknown reward probability of yielding a reward, where the objective is to pick whose machine's arm to play in order to maximize the accumulated reward over time. The dynamic multi-armed bandit (DMAB) [26] is the version of the MAB for dynamic scenarios, where the arm's reward probabilities can drift.
In the DMAB+MA, described in Algorithm 2, the DMAB was implemented as an adaptive operator selection approach to determine the best configuration of operators for a MA, that otherwise, follows the structure previously described in Section IV-A. The pool of MA operators [25] consists of four selection operators (stochastic, roulette, random and binary tournament), two operators for crossover (cyclic and orderbased), three mutation operators (insertion, reduced 3-swap and cumulative swap) and two survival strategies ((µ, λ) and (µ + λ)). Additionally, there are two fitness evaluation functions: the conventional function depicted in (1) and an alternative fitness function f , designed to deal with the intrinsic high neutrality of the CBSP which was previously reported in [27]. Each operator configuration is equivalent to a DMAB arm a i .
The DMAB+MA has been proven to be very effective, when compared with the MA and other existing reference algorithms from the CBSP literature [25]. The DMAB begins by playing all arms and assigning them initial reward values. It employs the upper confidence bound (UCB1) [62] to assign to each arm a i a confidence estimation based on its historical rewards and the number of times they have been used. This allows to balance the choice between arms that produced good results and those that have been underused. At each iteration the arm with the highest confidence is played and a reward on function of the fitness improvement it was able to achieve is assigned to it. The Page-Hinkley test (PH-Test) [63] is implemented to help the DMAB algorithms to adjust to scenarios were the most suitable arm to play can change over the time. When this occurs, the PH-Test triggers and all the existing arms are restarted.

A. EXPERIMENTAL SETTINGS
A set of 23 well-studied CBSP instances, listed in Table  2, was selected. The set contains topologically diverse, and representative graphs of order 24 ≤ n ≤ 200, with a number of edges ranging from 68 to 2000. The topologies included in this set are: 2 paths, 2 cycles, 2 wheels, 6 Cartesian products, 4 power of cycles, and 7 Harwell-Boeing graphs.
For constructing the STNs, we record the transitions between representative solutions that occur when the MA and DMAB+MA algorithms, described in Section IV, are used for solving the selected CBSP instances. To this end, MA and DMAB+MA were coded in C language and compiled in gcc 4.4.7 using an Intel Xeon CPU X5550 at 2.67 GHz Algorithm 2: DMAB Algorithm 1: input Guest graph G 2: output Best-found solution g 3: P ← initializePopulation(P ) 4: P best ← best individual in P 5: g ← P best 6: Set confidence and number of times arms have been played to zero 7: for i ← 1 to k do 8: P ← playArm(a i , P ) 9: Assign initial reward to a i 10: end for 11: repeat 12: Compute confidence for all arms 13: a s ← selectArm() 14: P ← playArm(a s , P )

B. MULTIDIMENSIONAL SCALING VISUALIZATION
After obtaining the samples from the MA and DMAB+MA algorithms, the pair-wise interchange distances [54] between unique solutions are computed and recorded in a distance matrix D. Then, our proposed partitioning method, based on multidimensional scaling [28], [29], maps the sampled solutions (in n dimensions) to 3-tuples in the Euclidean 3D space represented with an n × 3 matrix, called X . This is done by employing double centering and eigenvalue decomposition over the matrix D. The resulting coordinate matrix configuration X minimizes a loss function, called strain, which warranties to find a mapping respecting as much as possible (i.e., with the smallest error rate) the inter-location distances stored in D. All along this procedure, we keep track of the algorithm that produced each solution. Figures 1 and 2 present the results obtained by our modeling and partitioning steps, based on multidimensional scaling, when applied to two representative instances: can144 and path200 (similar results are obtained with the rest of the tested graphs). The fitness of each point in these plots corresponds to the average of the fitness values for all the solutions that share the same Euclidean position. A color map is used to represent with dark violet, solutions with small (better) fitness values, and with light yellow tones those having large fitness costs. This kind of visualization, previously used in [30], [55], is useful to reveal the fitness variations around specific known solutions, such as the bestfound. However, they lack the trajectory component to reveal how the transitions among solutions occurred through the search process.
For instance, in Figure 2 we clearly observe at right hand different zones of the search space which concentrate reduced  groups of solutions with small fitness values (dark violet points), while at left hand a cluster of yellow points (i.e., large fitness cost points) is surrounded by medium fitness cost points. The former represent a zone of the search space less frequently visited by the algorithm, while the latter one seems to be a more accessible region for the algorithms. From these plots, nonetheless, it is not possible to identify if the analyzed algorithm has followed a particular search trajectory connecting these two distant zones of the space.
As we will see the STNs, presented next, represent a better alternative visualization that overcomes this drawback.

C. STN ANALYSIS
The following step in constructing the STNs consist in recovering the set of search space sites visited by each algorithm (i.e., the nodes of the STN) from the coordinate matrix X . Next, the obtained coordinate values are rounded up to integers, and those sites having exactly the same resulting 3D coordinates are grouped to become a location of the STN under construction. At that point, these locations are used to produce two single STNs: one for each analyzed algorithm. Both single STNs are combined to produce a merged STN ST N Ai,Aj = (V Ai ∪ V Aj , E Ai ∪ E Aj ), which allow us to compare the search dynamics of the two best-known algorithms for solving the CBSP.
The set of steps described previously, to produce merged STNs, were applied to each one of the 23 selected CBSP instances in this experiment. Figures 3 and 4 present the resulting STNs for four representative host graphs: ibm32, can144, p9c9, and path200. In these figures we can observe at right hand merged STNs produced by employing our proposed partitioning method. In contrast, those STNs depicted at left hand were built up without mapping single solutions to locations in the search space, i.e., without applying our multidimensional scaling based partitioning method. The subnetworks colored in blue correspond to the MA algorithm, while those in orange belong to the DMAB+MA algorithm. The first nodes in the search trajectories are marked in yellow, the best ones in red, and those visited by both analyzed algorithms (shared nodes) in light gray. The size of the nodes is proportional to the number of times it was visited. Additionally, the end nodes with worse fitness value than the one of optimal/best-known solutions are depicted in dark gray color to identify them. Table 3 presents the set of test instances we employed in this experimentation, including the name, order, size, and value of the optimum/best-known solution [24], [25] for each one of them. It continues by presenting network metrics evaluated for the single STNs of each algorithm analyzed (MA, and DMAB+MA), regarding their number of nodes, length of the shortest path to the closest end node/best node, and the number of paths of such length. Shortest paths of length zero correspond to cases where the initial solution and the optimal/best-know solution were close enough to be mapped into the same network location by our search space partitioning method. The structural metrics for the merged STNs are introduced in Table 4. It displays the number of nodes and edges, the number of shared nodes, best nodes, and end nodes. The in-strength for end nodes is marked in bold when it is equal or larger than the in-strength of the best nodes. Problem instances where this occurs can be considered quite hard to solve, since it indicates that the search trajectories can be deviated to locations of inferior fitness quality.
One of the first noticeable aspects in the STNs produced is that the trajectories of the DMAB+MA algorithm are shorter and less likely to visit previously encountered solutions (to produce search cycles). They finish more frequently in best nodes and the DMAB+MA subnetworks, in the merged STNs, have fewer nodes as well as shorter paths. This supports previous claims for DMAB+MA to be the best from the two studied metaheuristics, as it clearly demonstrates that the algorithm can obtain better results, in terms of solution quality, than MA while offering more efficient search space exploring capabilities.
Our STN analysis also shows that the trajectories for both analyzed algorithms tend to finish in very different regions of the search space. For example, in the STN constructed for instance can144, shown in Figure 3(c), there are eleven unique locations corresponding to optimal/best-known solutions (red nodes), and nine other that are end nodes of worse quality (dark gray nodes). In general, the end and best nodes tend to not be shared and to belong to distinct regions of the search space, even after applying the partitioning process as it can be seen in Figure 3(d). This is an indicator of a multimodal fitness landscape with high quality solutions scattered across it.
The shared nodes in the STNs provide interesting information, since they represent regions explored by both algorithms and the ways on which they diverted from there. A high number of shared nodes could be associated to the existence of narrow funnel structures in the underlying fitness landscape. In Figures 3(b), 3(d), 4(b), and 4(d) it was observed that shared nodes tend to be close to the best and end locations. In several cases, only the DMAB+MA algorithm is able to reach the best nodes after passing from a shared node. Furthermore, DMAB+MA is able to accomplish this task employing fewer transitions.
The problem instances constructed as the Cartesian product of two graphs (p9p9, c9c9, k9k9, p9c9, p9k9, and c9k9) usually have a higher number of shared nodes in proportion to the total number of nodes in the corresponding STN (compare columns 2 and 4 in Table 4). This could signal that the fitness landscapes of these problem instances have structures that conduct to different high quality areas.
The results of the STN analysis using graph metrics, reported in Tables 3 and 4, reveal interesting information related to the structure of the networks and the differences in difficulty of the graph topologies tested. The instances of the path, cycle, and power of cycle topologies had a single best node regardless of their size (see column 5 in Table 4). Since the value of the optimum is known for these topologies, we know that the network locations corresponding to those best nodes contain in fact optimal solutions.
In the path and cycle topologies, the in-strength of the best node becomes smaller than that of the end nodes for the larger instances (see columns 7 and 8 in Table 4). This suggests that the search space for the path and cycle graphs contains very few optimal solutions. Finding the optimum then becomes harder because it seems to be also larger areas of suboptimal VOLUME 4, 2016 path100  99  304  11  9  376  6  9  path200  199  378  19  10  214  11  8  cycle100  100  343  3  10  363  15  10  cycle200  200  260  8  10  218  8  8  wheel100  2600  60  2  31  383  1  24  wheel200  10200  101  5  21  223  2  11  cyclePow100-2  300  266  2  10  303  2    solutions that result difficult to escape from, specially for the MA. In this sense, these two topologies could be considered harder to solve than the power of cycle one, since the instrength for the best nodes was always higher than for the end nodes. In fact, their in-strength of the best node increased for those instances with power k = 10, which have a larger number of edges. The wheel graphs had smaller in-strength for the best nodes when the instance size increases as well, however their number of different best nodes is larger. They present also more paths leading to them, which would imply that this topology is easier to tackle, because the optimal solutions are not as isolated as in the case of the cycle and path graphs.
The group of the Cartesian products (p9p9, c9c9, k9k9, etc.) presented some of the highest numbers of best nodes, generally in combination with numerous shortest paths to-wards them and a small number of end nodes. While this was also the case for some of the Harwell-Boeing graphs, like can24 and nos4, this subset of instances is not as homogeneous as the other ones, since it contains problems from diverse engineering areas. Therefore, we observed more variation in the structures of their STNs, for example, curtis54 and will57 both had very small in-strength for the best nodes and for can144 it was equal to that of the end nodes.

VI. CONCLUSIONS
This work presented an analysis, based on search trajectory networks (STNs), of two sate-of-the-art metaheuristics for the cyclic bandwidth sum minimization problem (CBSP): MA [24] and DMAB+MA [25].
Search space partitioning is an essential step during the construction of STNs models, which consists in mapping (c) STN for instance path200, without partitioning.
(d) STN for instance path200, after partitioning. the solution sampled by the studied algorithms to locations in the search space. We introduced a novel search space partitioning method based on the classical multidimensional scaling [28], [29] for reducing solutions (in n dimensions) to 3-tuples in the Euclidean 3D space, while preserving their distances. This partitioning method is more generally applicable than the previously proposed methods [15], [16], as it can be applied to any solution representation, either discrete or continuous, as soon as a suitable distance function between solutions is available.
The STN analysis carried out helped us to infer different important characteristics of the fitness landscape associated to the problem of study. For some of the CBSP instances we considered (wheel, Cartesian product, and Harwell-Boeing graphs), their search space contains multiple optimal/bestknown solutions that are sparsely distributed across the fitness landscape. Meanwhile, for other graph topologies (path, cycle, powers of cycle), there seem to be specific regions containing the optimum. It also allowed us to identify that certain regions traversed by both studied algorithms are close to high quality areas, but it is usual that only the DMAB+MA algorithm has the ability to reach them, while the MA instead gets to worse quality end nodes.
The evidence from this STN-based study helps to demonstrate that DMAB+MA is more efficient for conducting the search process, as observed in its shorter trajectories, avoidance of areas where the MA gets trapped, shorter paths to the optimal/best-known, and higher success rates.
Further work would include: a) refining and further testing the partitioning method for other combinatorial problems, b) using STNs to analyze how the change in the fitness function [27] component alone helps to shape the search trajectories, and c) researching the correlation (if any) between the STN measurements and the actual size of the attraction basins.
GABRIELA OCHOA is a Professor of Computing Science at the University of Stirling in Scotland. Her research lies in the foundations and applications of evolutionary algorithms and metaheuristics, with emphasis on autonomous search, fitness landscape analysis and combinatorial optimisation. She holds a PhD from the University of Sussex, UK, and has held academic and research positions at the University Simon Bolivar, Venezuela, and the University of Nottingham, UK. Her recent work on network-based models of fitness landscapes has enhanced their descriptive and visualisation capabilities, producing a number of publications including 4 best-paper awards and 4 other nominations at leading venues. She collaborates cross-disciplines in the use of computational intelligence in healthcare and conservation. Her Google scholar citations currently report over 5,600 citations with an H-index of 38. She has been active in organisation and editorial roles within leading Evolutionary Computation venues including the IEEE Transactions on Evolutionary Computation and the Evolutionary Computation Journal. Since 2008, he has been an Associate Professor (Cinvestav-3B Researcher) with the Cinvestav -Tamaulipas, Victoria, Mexico. He has authored and co-authored a book and over 50 technical papers and book chapters. His publications currently report over 690 citations in Google Scholar with an H-index of 14. His current research interests include evolutionary computation as well as the design and implementation of effective metaheuristic algorithms for solving large-scale combinatorial optimization problems arising in various application areas like bioinformatics and graph theory.