1 Introduction↩︎

The Traveling Salesperson Problem (TSP) is a classical \(\mathcal{NP}\)-hard optimization problem of utmost relevance, e.g., in transportation logistics, bioinformatics or circuit board fabrication. The goal is to route a salesperson through a set of cities such that each city is visited exactly once and the tour is of minimal length. In the past decades tremendous progress has been made in the development of high-performing heuristic TSP solvers. The local search-based Lin-Kernigham Heuristic (LKH) [1] and the genetic algorithm Edge-Assembly-Crossover (EAX) [2], along with their respective restart versions introduced in Kotthoff et al. [3], undeniably pose the state-of-the-art in inexact TSP solving.

Automated Algorithm Selection (AS), originally proposed by Rice [4] back in 1976, is a powerful framework to predict the best-performing solver(s) from a portfolio of candidate solvers by means of machine learning. It has been successfully applied to a wide spectrum of challenging optimization problems in both the combinatorial [5]–[9] and continuous domain [10], [11] with partly astonishing performance gains – see the recent survey by Kerschke et al. [12] for a comprehensive overview. In particular, the TSP was subject to several successful AS-studies [3], [13]–[16] which exploited the complementary performance profiles of simple heuristics on the one hand and the state-of-the-art solvers LKH and EAX on classical TSP benchmark sets on the other hand.

In the classic setting, AS relies on characteristic problem instance properties, termed (instance) features. These features are used as predictor variables for classical machine learning algorithms, e.g., random forests or support vector machines. The key idea – and ideal outcome – is that these features can easily be used to automatically derive decision rules that are well-suited to partition the instance space into ideally disjoint sub-spaces of instances, and which then are uniquely solved best by different solvers. However, features have many drawbacks: they are usually hand-crafted in a tedious process [17], partly require time-consuming calculations (which need to be taken into account by the model fitting step) and are problem-tailored (or at least specific to a problem domain, e.g., graph problems). Moreover, we usually prefer light-weight models with few features. Hence, training AS models is frequently combined with automated feature selection methods [18]–[20] or dimensionality reduction techniques [21].

Recently, Alissa et al. [22] took a promising new path after exploring first approaches to avoid manual feature-crafting [23], [24]. The authors proposed a deep learning based approach which does not rely on any a-priori calculated feature profiles. Instead, in their study on the 1D-Bin-Packing Problem (BPP) the neural network is given temporal sequence data of the BPP instance as the only input. They were able to achieve drastic improvements. In this paper we adopt and adapt this idea for the TSP. To this end, we evolve a set of instances where LKH and EAX show strongly different behaviour in terms of Penalized Average Runtime¹ (PAR10; [25]). We show that with classical AS we can clearly beat the Single Best Solver (SBS; the solver with best average performance). However, the gap to the Virtual Best Solver (VBS; perfectly predicting oracle) can only be reduced slightly with much room for improvement. This holds true even in the case when we enrich the machine learning pipeline with (1) hand-selected feature subsets (based on exploratory data analysis), (2) different feature selection methods, or (3) a combination of both. After that, we propose a feature-free deep learning approach where the neural networks are trained on the plain image representations of Euclidean TSP instances. This approach achieves competitive performance, but drops the need for manual feature derivation and calculation.

The remainder of this paper is structured as follows. We describe the benchmark set (generation) and pre-selection of feature subsets in Sections 2 and 3 respectively. In Section 4 we present the results that we achieve using classical feature-based AS approaches. Next, in Section 5, we detail our feature-free deep learning based approaches, and compare the results with the classical models. We close with a discussion and an outlook in Section 6.

2 Evolving TSP Instances↩︎

Our benchmark requires a set of Euclidean TSP instances that show strong differences in algorithmic performance. To this end, we adopt an evolutionary algorithm (EA) and creative mutation operators recently proposed by Bossek et al. [26]. Their method allows for the tailored generation² of TSP instances that (1) have the desired performance difference (in terms of the ratio of PAR10-scores), (2) show multifarious topologies in terms of point arrangement in the Euclidean plane³, and (3) are well distributed in the space of instances characteristics/features; in particular the latter two properties were not achieved by evolutionary instance generation methods before [14], [15]. For sake of brevity we refer the reader to Bossek et al. [26] for more details.

We generated a balanced data set of \(1\,000\) TSP instances with \(n=1\,000\) nodes per instance using the EA parameters from [26]; each \(500\) being uniquely faster solved to optimality by either EAX or LKH.⁴ All generated data is available in a public GitHub repository (https://github.com/mvseiler/PPSN_TSP_DL). Fig. 1 depicts the performance – measured by means of PAR10 – of EAX and LKH on the entire benchmark set. The plot highlights apparent potential for automated algorithm selection due to strong performance differences. On the other hand the data reveals the general superiority of the EAX solver since it is much harder to evolve instances that are hard for EAX (shown as orange points). There are just two instances for which the EAX hits the cutoff time \(T=3\,600\) seconds (1 hour) at least once in all of its ten independent runs on that instance. In contrast, LKH frequently gets stuck in local optima – see the cluster of green points between \(T=3\,600\) and \(10 \cdot T\) in the top left corner.

3 Identifying Adequate Subsets of TSP Features↩︎

The state-of-the-art TSP-related feature sets [14], [16], [28] consist of hundreds of hand-crafted features. Features range from statistics (mean, variance etc.) of edge lengths, angles of nearest neighbors, to more sophisticated features based on Minimum Spanning Trees (MST) or \(k\)-Nearest-Neighbor-Graphs (\(k\)-NNG). Fig. 2 depicts visual impressions of MSTs and \(k\)-NNGs on two evolved instances; these images will be a key ingredient to the neural network in Section 5.

Due to the size of feature sets, in the context of algorithm selection, automated feature selection methods have shown their suitability for automatically filtering a (small) subset of discriminating features [13]. Regardless of the sophistication of the feature selection method at hand, feature selection needs to cope with an exponentially sized search space of the underlying subset-selection problem. Hence, in order to assist our feature-based model fitting we conduct a simple univariate exploratory data analysis in order to identify a lucid subset of adequate features a-priori. To this end we adopt a simple greedy heuristic. First, all features \(f\) are scaled to \([0,1]\) to allow for a fair comparison across the features. Then, we perform a two-sided non-parametric Wilcoxon-Mann-Whitney test [29] at significance level \(\alpha=0.05\) per feature \(f\), to check the null hypothesis that the distributions of the EAX instances and the LKH instances with respect to \(f\) differ by a location shift equal to zero. We extract the 15 most relevant features according to the smallest \(p\)-values. This procedure was done once for the entire benchmark set of salesperson features as the most comprehensive feature set (see Sec. 4) and was repeated for subsets of “hardest” instances of decreasing size. Hardest in this context relates to (a) each 300 and 150 hardest instances for each solver with respect to mean PAR10-score, or (b) the ratio of mean PAR10-scores. Both methods pursue to reduce the benchmark set to instances with maximal performance differences.

Fig. 3 shows the distribution of selected features both across the whole benchmark set and the subset of each \(300\) most difficult instances with respect to mean PAR10-score for each solver (marked with crosses in Fig. 1). Noticeably, the 15 most relevant features are identical for both sets. They are mainly composed of summary statistics on strong connected components of the nearest neighbor graph (nng_*) and properties based on minimum spanning trees (mst_*). This is very much in line with crucial features identified in most TSP-related AS-studies [13]–[16], [30], [31] by sophisticated variable importance measurement. These features seem plausible since both MSTs and NNGs capture the global structure, e.g., existence of clusters etc., very well. A close look at Fig. 3 suggests, that instances that are easy for EAX cover a wider range of values which is derived from wider (green) boxes for the features while easy instances for LKH show much more narrow (orange) boxes, with many strong outliers though. The right hand plot, however, shows that the hardest instances seem to be better separable with the features. For instance, for the features in the top 4 rows we observe that \(75\%\) of feature values for LKH-friendly instances are higher than \(75\%\) of the respective values for EAX-friendly instances.

4 Classical Algorithm Selection↩︎

Fig. 1 reveals very complementary performances of EAX and LKH, which gives us reason to assume that automated AS might work well in this setting. Further, as outlined in the univariate analysis of the TSP features (see Sec. 3), the features also indicate their potential for distinguishing instances that are beneficial for EAX from instances for which LKH is preferable. As previous works [3], [13] already confirmed the effectivity of feature-based AS, we adopted their experimental setup – and only slightly modified it to the scenario at hand. Below, we will outline the considered machine learning algorithms, feature sets, as well as feature selection strategies, which have been used for training our final selectors.

4.1 Experimental Setup↩︎

All our candidate AS models are trained using PAR10 [32] as performance measure and assessed with a 10-fold cross-validation (CV). As indicated by Fig. 1, it is much more difficult for LKH to perform well on the instances that were evolved in favor of EAX (green points), rather than vice versa (orange). Therefore, the selectors will likely have a bias towards EAX instances. To adjust for this bias, we additionally tune the classification threshold for all trained models. For training the potential automated AS models, we considered four different classifiers [33] using the R-package mlr [34]: decision trees [35], random forests [36], support vector machines [37] and gradient boosting [38]. Each of them is trained using three different feature sets: the UBC features from Hutter et al. [28], a large set of features by Pihera and Musliu [16], as well as the TSP features from the R-package salesperson⁵ [39]. The salesperson features provide the up-to-now most comprehensive collection of features; in fact, they are a strict superset of the Pihera and tspmeta features [15]. On the other hand, the UBC and Pihera features led to the best performing algorithm selectors in previous works [3], [13] – which did not consider the salesperson features as the package did not exist back then. Note that there is a large overlap across the three considered feature sets as outlined in [12]. To reduce the noise within and redundancy between the features, we additionally created five small subgroups from the salesperson feature set, consisting of 15 features each (see Sec. 3 for details).

In addition to the 32 potential AS models described above (8 feature sets \(\times\) 4 learners), the respective feature sets were further reduced using three automated feature selection strategies: sequential floating forward selection (sffs), sequential floating backward selection (sfbs), and – for the reduced feature sets – exhaustive search of the 15 features [12], [40]. This resulted in 84 further candidate selectors.

4.2 Findings↩︎

Tab. 1 summarizes the averaged PAR10 performances of all 116 considered AS models, with the best achieved scores highlighted in red. Of course, all shown PAR10-scores already include the costs for the computation of the TSP features. On average those costs account for merely 0.7s at most. According to the listed performances of the best models (61.21s), the SBS (67.47s) and the VBS (4.92s), the best classical AS approaches are able to reduce the SBS-VBS-gap by 10%.

The best found selectors are random forests, which reduced the top 15 features that they were given initially, to the following four features: the sum, arithmetic mean, median and coefficient of variation of the distances of the MST. Moreover, the tuned thresholds varied from 4% to 33% across the ten folds, implying that EAX has always been selected once the model predicted EAX with a probability of at least 33%. In consequence, out of all instances, in which LKH was actually the faster solver, it has only been selected 151 times (corresonding to roughly 30%). On the other hand, LKH has only been (wrongly) picked in 4% of the cases, in which EAX would have been the correct choice.

When starting with the full feature sets from Pihera, UBC and salesperson, a support vector machine based on a subset of the UBC features achieved the best performance (printed in bold in Tab. 1). In fact, the PAR10-score of 56.67s is only slightly worse than the one of our best selector(s). Noticeably, the SVM also relied on MST features only: the arithmetic mean and standard deviation of the lengths of the edges in the MST, as well as the skewness of its node degrees.

Our findings are also confirmed by the left image of Fig. 4 as only few observations (20) are located above the diagonal. However, looking at the right image, it becomes nearly obvious that the selector is still quite far away from the performance of the VBS, as shown by the many misclassifications (381 out of all 1 000 instances) above the diagonal.

In an attempt to better understand the reason for these rather small thresholds, we investigated the misclassification costs in detail. Prior to tuning the threshold, LKH was predicted 154 times when EAX would have been the correct choice – and each of those misclassifications caused (on average) an overhead of 4 423.89s. In contrast, the 128 cases in which EAX was predicted instead of LKH only came with an average penalty of 95.95s. After tuning the thresholds, each of the 20 wrong predictions of LKH caused only 290.04s – compared to the average penalty of 126.07s for the 361 wrong predictions of EAX. Thus, by being rather conservative and only predicting LKH in cases, where the model is highly certain, the selector was able to reduce the misclassification costs significantly.

As our results indicate, the common feature-based approaches are able to improve over the SBS. However, it is also noticeable that the currently available features still have a very hard time in extracting sufficient information from the TSP instance to reliably predict the better solver. Therefore, we will test the suitability of deep learning neural networks as an alternative or supporting means for automated algorithm selection.

5 Deep Learning Based Approach↩︎

As demonstrated in the previous section, feature-based AS methods can outperform the SBS. However, these models come with three major drawbacks: they (1) are hand-crafted in a tedious process, (2) partly require time-consuming calculations, and (3) are problem tailored (see Sec. 1). To overcome these issues, we propose a novel, and sophisticated feature-free approach that is based on so-called Convolutional Neural Networks (CNN) [41].

5.1 Experimental Setup↩︎

To train CNN based AS models that are independent of the commonly used TSP features, we will produce different visual representations of our TSP instances (see Fig. 2) and use them for training the deep learning networks. Those images are created with a resolution of \(512 \times 512\) pixels and the coordinates of the instances are scaled to fill out the entire image as exemplarily shown in the two images in the left column of Fig. 2. In addition to these point clouds, images of corresponding Minimum Spanning Trees (MST, second column of Fig. 2) and \(k\)-Nearest-Neighbor Graphs (\(k\)-NNG, right column of Fig. 2) with \(k = 5\) were generated. We chose MST and \(5\)-NNG as additional visual representations because we found in Section 3 that the 15 most important features are based almost exclusively on MST and \(5\)-NNG graphs. In the \(5\)-NNGs, not only mutual (strong) connections but also one-sided (weak) links were considered, in which one city belongs to the nearest neighbor set of another city, but not vice versa.

Admittedly, only networks whose generation was based exclusively on point clouds can be described as feature-free. For a better comparison, however, we have additionally evaluated feature-based networks that were trained with images of the corresponding MST and \(5\)-NNG. Hence, we considered two different scenarios for our network-based approaches. In the first scenario (S1), the networks were trained based on (a) point clouds of the cities (Points), (b) MST images, and (c) \(5\)-NNG images. In scenario (S2), we combined (a) the scatterplots with the MST images (as two input channels), and (b) the scatterplots with the MST and \(5\)-NNG images (as three input channels). As the costs for generating the images are insignificant, we have not taken their generation time into account when computing the PAR10 scores of the deep learning models. For larger instances, though, these times would have to be taken into account.

To process the visual representations of the instances, we used eight stacked convolutional layers (see Fig. 5). Three of them used \(Strides=2\) to reduce the size of the feature maps and four of them used \(Dilation=\{2,3\}\)-Kernels to enlarge the receptive fields and, thus, gain a larger view of the instances (see Fig. [fig:cnn95working] to compare the effects of Strides and Dilation). We used Rectified Linear Unit (ReLU) [42] as activation function for all layers except for the last linear layer, for which we used a Softmax activation. To improve the training speed, the outputs of all convolutional layers are normalized by using Group Normalization (GN) [43] with \(G=8\). The GN layers are in-between the convolutional layers and the ReLU activation. For transition from the three-dimensional convolutional (\(Width \times Height \times Channels\)) layers to the one-dimensional linear layer, a Global Average Pooling Layer (GPL) [44] is used. After the GPL, a Dropout (DP) [45] layer with 25% dropout is added to improve regularization. The final layer is a single, linear layer with two output neurons – one for EAX and one for LKH (see Fig. 5). Last, we used 10-fold cross-validation to evaluate the performance of the neural networks. The folds were the same as for the classical feature-based approach. All networks were trained using mini-batches of eight, Adam [46] as optimizer and Cross-Entropy [47] as loss function. Note that neural networks are most commonly trained using Stochastic Gradient Descent [48], which strongly differs from the training methods used in Section 4.

5.2 Findings↩︎

The best performing classical approach achieved a mean PAR10-score of \(56.29\) (see Tab. 1). Our feature-free networks, which were trained exclusively using the points, achieved a mean PAR10-score of \(56.31\) after tuning the threshold and thus a similar performance (see Tab. 2) as the feature-based, classical approaches.

As stated before, we additionally investigated whether adding additional features to the networks could improve the models’ performances. Therefore, we also trained feature-based networks using MST and 5-NNG images. As shown in Tab. 2, both variants perform noticeably better. Besides, we found that while the thresholds between the points and the MST models are rather similar, the thresholds of the NNG models are on average \(10\%\) higher. Next, the thresholds of the Points and MST models range from \(13\%\) to \(34\%\), and \(16\%\) to \(35\%\), respectively, while the thresholds of the NNG models range from \(3\%\) to \(73\%\). Thus, networks trained on the NNG images appear to be less stable.

Moreover, the networks based on the Points correctly predicted EAX in 91.8% (and thus 459 times) of the cases, in which EAX was the better solver, compared to only 22.6% (113 cases) for the LKH-friendly cases. This behavior likely results from the fact that a misclassification of an instance, which is favorable for LKH, is cheaper than a misclassification of an instance that is easier for EAX. In contrast, the MST-based networks predict EAX in 91.6% (458) and LKH in 27.2% (136) cases, correctly. Thus, compared to the feature-free networks, which were exclusively based on Points, the MST networks benefit from correctly identifying LKH-friendly instances – without losing accuracy on the EAX-friendly instances. Noticeably, in case of the 5-NNG networks, only 84.6% (423) of the EAX-easy instances are classified correctly, compared to 34.8% (174) among the LKH-easy instances. Thus, despite the improvements among the instances that are favorable for LKH, the PAR10 score of the NNG-based networks is inferior to the MST-based selector, as misclassifying EAX-easy instances is more expensive.

To investigate whether the combination of the three different input variants would lead to networks that achieve better performances when predicting EAX- and LKH-easy instances, we combined Points and MST, as well as Points, MST and NNG into two and three input channels (Scenario 2), respectively. However, as shown in Tab. 2, combining the visual representations does not improve the networks’ overall performances. Also, while the network based on Points + MST classifies 31% (155) of the LKH-easy instances correctly, the selector based on Points + MST + NNG only succeeds in 17% (85) of the respective cases. We further observed that the threshold values are quite similar to the ranges of the Points and MST models from scenario S1.

As visualized in Fig. 6, the Points (left) and MST networks (right) prefer EAX over LKH – as there are far more observations below the diagonal. This is also confirmed by the low thresholds (see Tab. 2). Interestingly, the networks solely based on visual representations of the instance, even perform slightly better than the classical feature-based AS models – but are still clearly inferior to the VBS.

6 Conclusions and Outlook↩︎

The conducted experiments shed light on still existing shortcomings of classical feature-based per-instance algorithm selection on the TSP. While previous studies clearly reported successful approaches, the informative character of existing TSP feature sets reaches its limits for instances specifically evolved for maximum performance difference of the two state-of-the art heuristic solvers. Sophisticated mutation operators here lead to so far unobserved topological structures. Despite outperforming the SBS, the gap to the performance of the oracle-like VBS cannot be closed substantially, even after utilization of sophisticated preprocessing and feature selection approaches.

However, it again becomes obvious that the minimum spanning tree and nearest neighbor structures of the points are most informative in discriminating solver performances. We build on this information and enrich a deep neural network approach based on images of the instance topology by specific images visualizing the minimum spanning tree and the nearest neighbor graph. Most interestingly, our feature-free deep neural network nicely matches the performance of the quite complex classical AS approach (see Tab. 3), despite being solely based on an image of the instance’s points.

This proof-of-concept study thus shows the huge potential of deep learning, feature-free approaches in this domain which we will exploit in future studies by more sophisticated networks and altered loss functions specifically adapted to the designated performance indicators. Moreover, additional image channels will be added, e.g., in terms of heatmaps. Specific investigations have to be conducted with regard to the scaling behaviour of the approach, as image resolutions most probably will have to be carefully adapted to increasing instance sizes.

On the other hand, the observed limitations of classical TSP features show the necessity of enriching the library of TSP features by alternative sets which capture other kinds of — obviously important — instance structures. We will apply our approach to classical feature sets such as RUE or TSPLib as well for a comparison. However, it is specifically noteworthy that, in principle, the deep learning approach nicely generalizes to other graph-based optimization problems while instance features are almost exclusively tailored to the focused domain.

Acknowledgments↩︎

The authors acknowledge support by the European Research Center for Information Systems (ERCIS).

References↩︎