{'fraction_non_alphanumeric': 0.030673750948787784, 'fraction_numerical': 0.013707192927624236, 'mean_word_length': 4.492949110974862, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generalization of Machine Learning for Problem Reduction: A Case Study on Travelling Salesman Problems

Combinatorial optimization plays an important role in real-world problem solving. In the big data era, the dimensionality of a combinatorial optimization problem is usually very large, which poses a significant challenge to existing solution methods. In this paper, we examine the generalization capability of a machine learning model for problem reduction on the classic traveling salesman problems (TSP). We demonstrate that our method can greedily remove decision variables from an optimization problem that are predicted not to be part of an optimal solution. More specifically, we investigate our model's capability to generalize on test instances that have not been seen during the training phase. We consider three scenarios where training and test instances are different in terms of: 1) problem characteristics; 2) problem sizes; and 3) problem types. Our experiments show that this machine learning based technique can generalize reasonably well over a wide range of TSP test instances with different characteristics or sizes. While the accuracy of predicting unused variables naturally deteriorates the further an instance is from the training set, we observe that even when solving a different TSP problem variant than was used in the training, the machine learning model still makes useful predictions about which variables can be eliminated without significantly impacting solution quality.


Introduction
In the big data era, we are often confronted with optimization problems with thousands or even millions of decision variables, e.g., social network analysis (Balasundaram et al., 2011;Gao et al., 2018). The large problem size poses significant challenges to existing solution algorithms, especially to generic Mixed Integer Programming (MIP) solvers such as CPLEX, which typically has difficulty in optimally solving or even finding good solutions for such large-scale optimization problems in a reasonable computational time. Moreover in many practical applications, e.g., trip planning (Friggstad et al., 2018), we need to provide a high-quality solution to users within a few seconds. This is hard to achieve especially when the problem size is very large, which necessitates the use of an effective problem reduction technique that can significantly prune the search space but still capture an optimal (or near-optimal) solution in the reduced space.
Recently, there has been a growing trend of applying machine learning for problem reduction (Li et al., 2018;Grassia et al., 2019;Ding et al., 2019;Sun et al., 2019). These machine learning models are typically trained on easy problem instances for which the optimal solution is known, and predict for a given hard unsolved problem instance a subset of decision variables that most likely belong to an optimal solution. By greedily removing decision variables that are not expected to be part of an optimal solution, the original large search space can be significantly reduced to a size that is manageable by an existing solution algorithm. In our recent work (Sun et al., 2019), we have developed such a Machine Learning model for Problem Reduction (MLPR), which builds on statistical measures computed from stochastic sampling of feasible solutions. We have empirically shown that as a prepossessing technique, our MLPR method can significantly improve the performance of existing solution algorithms when used to solve large maximum weight clique problems.
Although the idea of problem reduction using machine learning is generic, it is still unclear whether our MLPR method is also effective on combinatorial problems other than maximum weight clique problems. In this paper, we examine the effectiveness of our MLPR method on the classic traveling salesman problem (TSP). We consider TSPs on a complete graph where the objective is to search for a shortest route that visits each vertex and returns to the original vertex in the graph. We adapt our MLPR model to predict for each edge whether it belongs to a shortest route, and remove from the complete graph those who do not. The aim is to find a sparse subgraph that still contains a (near) optimal tour. This adaptation is nontrivial, because problem-specific features and sampling methods have to be designed for TSPs. Furthermore we parallelize our MLPR model in this paper, so that the computational time of our MLPR model can be significantly reduced by using multi-cores.
Since the TSP has been extensively studied and many effective solution algorithms have been developed for solving TSPs (Applegate et al., 2006a,b;Lin & Kernighan, 1973;Helsgaun, 2000), our primal goal here is not to further push the limit of problem solving. Instead, we focus on exploring the generalization capability (more specifically generalization error) of our MLPR model when training and test instances are different. Generalization error is very relevant to real-world problem solving, because in practice the test instances on which a trained model is applied are potentially quite different from the training instances. For example, the routing problem that a navigation company solves on a regular basis might drift over time.
A significant contribution of this paper is to provide a systematic analysis on the robustness of our MLPR model when such nontrivial changes happen in test instances. We empirically show that our MLPR model generalizes reasonably well to a wide range of test TSP instances with different characteristics or sizes. We also identify where our MLPR model may not perform well, i.e., on the test instances that are deliberately made to be very different from the training instances in terms of problem characteristics. This provides guidance on how to construct a good training set and when to update the training set in practical contexts.
Taking a step further, we investigate whether the knowledge learned from one variant of TSP instances can be transferred to solving other TSP variants. Our experimental results show that the MLPR model trained on symmetric TSP instances performs fairly well on some of the test instances from other TSP variants, although we do observe a performance degradation when the test TSP variants are gradually moved away from the training TSPs. This indicates it is possible to develop a more generic MLPR model that does not require re-training when applied to different problems (or at least a class of problems).
The remainder of this paper is organized as follows. In Section 2, we briefly describe the background and methods related to problem reduction. In Section 3, we adapt our MLPR model to reduce problem size for TSP. Section 4 presents the experimental results. The last section concludes the paper and suggests potential research directions for future work.

Background and Related Work
We briefly describe TSP in Section 2.1, and review the problem reduction techniques based on machine learning in Section 2.2. Because our MLPR model uses support vector machine (SVM) as the classification algorithm, we will briefly describe SVM in Section 2.3.

Travelling Salesman Problem
Given n cities {v 1 , v 2 , · · · , v n } and pairwise distance between cities {c i, j | i, j = 1, · · · , n, i = j}, the objective of the TSP problem is to find the shortest route that visits each city and returns to the original city. We use u i to denote the visiting order of city i, and use a binary variable x i, j to denote whether city j is visited directly after city i. Without loss of generality, we set u 1 = 1 (route starts from city 1). The Miller-Tucker-Zemlin formulation of TSP can be written as The first two constraints, (2) and (3), ensure that each city is arrived at and departed from exactly once; and the constraint (4) eliminates subtours. More computationally efficient formulations exist, but this is sufficient for logical correctness. Note that of the n 2 x i j variables, exactly n(n − 1) must be zero in any feasible solution. Removing such variables that are not part of any optimal solution would give a smaller problem with the same optimum. The TSP has been intensively studied and many solution algorithms have been developed to solve this problem, e.g., the Concorde exact solver (Applegate et al., 2006a,b), the Lin-Kernighan heuristic method (Lin & Kernighan, 1973;Helsgaun, 2000) and the "backbone" based heuristics (Dong et al., 2009;Jäger et al., 2014). Recently, there has been a growing interest in using machine learning to automatically learn a solution algorithm to solve combinatorial optimization problems (Bengio et al., 2018). The learning-based methods for solving TSP include: Vinyals et al.

Problem Reduction Based on Machine Learning
Many combinatorial optimization problems contain a large number of decision variables, most of which are "irrelevant" to the optimal solution. For example in a symmetric TSP with n cities, the total number of binary variables is n(n − 1)/2, while a shortest route only uses n binary variables (for those the value is 1). The goal of problem reduction is to identify some of these irrelevant variables and removes them from the original problem, in the hope that the reduced problem can be solved more easily. However, identifying these irrelevant variables is a nontrivial task itself.
Most of the existing problem reduction methods in mathematical programming are exact approaches, which only remove decision variables that cannot be part of an optimal solution, based on mathematical reasoning and/or computation of an objective bound (Jonker & Volgenant, 1984;Hougardy & Schroeder, 2014). An exact approach guarantees that the reduced problem always contains an original optimal solution, but in many cases it is computationally expensive and/or is not equipped with means to significantly reduce the problem size.
Fortunately, many combinatorial optimization problems have a "backbone" structure; that is the optimal solution of a problem is likely to share some components with high-quality solutions (Kilby et al., 2005;Wu & Hao, 2015). This makes it possible to statistically quantify which decision variables or solution components are more likely to be part of an optimal solution from sample solutions. This heuristic approach, although it does not have an optimality guarantee, can usually remove a large number of irrelevant decision variables from a given problem instance (Fischer & Merz, 2007;Sun et al., 2019).
Our MLPR method originally proposed in (Sun et al., 2019) belongs to the class of heuristic reduction approaches. We use optimally-solved problem instances as training set, and apply machine learning to automatically learn a rule to separate the decision variables that belong to an optimal solution from those who do not (irrelevant variables). We extract computationally-cheap problem features as well as statistical measures computed from random samples of feasible solutions to characterize each decision variable. Based on these features, we predict for each decision variable a likelihood of whether it belongs to an optimal solution. Our MLPR method can be used as a preprocessing technique to remove decision variables that are not expected to be part of an optimal solution from an unseen test problem instance. We will describe our MLPR model in more detail in Section 3.
Closely related to our MLPR method, Lauri & Dutta (2019) also developed a machine learning model for problem reduction to list all maximum cliques in a graph. They only use features directly computed from graph data to characterize a vertex, and remove vertices that are predicted not to be part of a maximum clique. Building on this work, Grassia et al. (2019) developed a multi-stage pruning technique that can further reduce problem sizes for sparse graphs. They also investigated the effects of removing edges instead of vertices from a graph. The main difference between these methods and ours is that these methods do not use statistical features computed from stochastic samples of feasible solutions. It is worth pointing out that these statistical features are of vital importance to our MLPR model, which helps our model generalize well to test problem instances that are not seen during training. We will describe these statistical features in Section 3.1.2.
Apart from using problem reduction as a preprocessing step (i.e., removing irrelevant decision variables from a given problem instance a priori), there exist other smart uses of problem reduction techniques. He et al. (2014) learned a node pruning policy for branch-and-bound algorithms to heuristically cut off branches that are unlikely to generate a better primal solution. Li et al. (2018) estimated a likelihood for each decision variable of whether it belongs to an optimal solution, and used the estimated probabilities to guide a tree search algorithm. Ding et al. (2019) trained a graph convolutional network to predict solution values for binary variables, and used the predicted values to generate a global inequality constraint to prune the search space. These methods are typically designed for a particular type of solution algorithms. In contrast, our MLPR method is more generally applicable and can be used as a preprocessing technique for any existing solution algorithm.

Consider a binary classification task with m training instances
where f f f i is the feature vector and l i ∈ {−1, 1} is the class label of the i th training instance. A classification algorithm aims to find a decision boundary to separate the positive (label 1) and negative (label −1) training instances as well as possible.
We first assume the positive and negative training instances can be separated by a linear classifier (h), parameterized by (w w w, b): where sgn(w w w T The geometric margin of (w w w, b) with respective to a training set S is the smallest geometric margin to the individual training instances: The aim of SVM (Boser et al., 1992;Cortes & Vapnik, 1995) is to find a decision boundary, determined by (w w w, b), that maximizes the geometric margin γ: SVM is also known as an optimal margin classifier. Scaling w w w and b by any positive number does not change the decision boundary: w w w f f f + b = 0. Thus we can restrict the norm of w w w to be any positive number without changing the optimal decision boundary. In order to efficiently solve the optimization problem, the norm of w w w is usually set to 1/γ, i.e., ||w w w|| = 1/γ. Thus maximizing γ is equivalent to maximizing 1/||w w w||, which is also equivalent to minimizing 1 2 w w w T w w w. The optimization problem is then transferred to a quadratic programming with linear constraints, which can be solved efficiently: Let (w w w * , b * ) determine the optimal decision boundary. The The support vectors are the training instances with the smallest geometric margin (those closest to the optimal decision boundary). Thus only support vectors can influence the optimal decision boundary. Adding or deleting a training instance which is not a support vector does not change the optimal decision boundary. When the training set S can not be well separated by a linear classifier, we can map the feature vector f f f to a higher-dimensional space using a non-linear function φ (·), in the hope that the training instances can be separated more easily in the higher-dimensional space. We can also use regularization, that allows a smaller geometric margin at a cost of increasing the objective value. Importantly this also caters for the case where the given set of non-linear functions is unable to provide a correct classification for all training instances. The primal optimization problem becomes where φ ( f f f i ) maps the feature vector f f f i into a higher-dimensional space; r + > 0 and r − > 0 are the regularization parameters for positive and negative training instances; ξ i , i = 1, · · · , m are slack variables and g(·) is a loss function. We will denote SVM with first order loss function g(ξ i ) := ξ i as L1-SVM, and SVM with second order loss function g(ξ i ) := (ξ i ) 2 as L2-SVM. If the function φ ( f f f i ) maps the feature vector f f f i to a very high dimensional space, solving the primal optimization problem is computationally slow. In this case, the dual problem may be easier to solve. Considering L1-SVM, the dual optimization problem is where {α 1 , · · · , α m } are dual variables, e e e = [1, · · · , 1] T is the vector of all ones, Q is an m × m positive semidefinite matrix, and is the kernel function. The kernel function avoids the need to explicitly compute φ (·), thus is computationally efficient. For example the radial basis function (RBF), defined as where γ k is a kernel parameter, implicitly maps the feature space to an infinity dimensional space. But computing the RBF kernel only the optimal decision boundary and thus the prediction for a given new instance is only determined by the support vectors (those with α * i > 0). It is noteworthy that the number of support vectors is usually much smaller than the number of training instances in the training set.

Problem Reduction for Travelling Salesman Problem Using Machine Learning
In this section, we adapt our MLPR method originally proposed in (Sun et al., 2019) to prune the search space for TSP. We model TSP as a complete graph G(V, E,C), where V denotes a set of cities, E denotes edges between cities, and C denotes edge costs (e.g., distance between cities). The objective of TSP is to search for a route with minimum edge costs that visits each vertex and return to the original vertex. Our MLPR method uses machine learning to predict for each edge whether it belongs to an optimal route, and removes from the complete graph the edges that are not expected to be part of an optimal route.
We use optimally-solved TSP instances as training set, and treat each edge in a solved graph as a training instance. We assign a class label 1 to the edges that belong to the optimal route and −1 to those who do not. We will extract two statistical measures and four graph features to characterize each edge in Section 3.1. After constructing a training set, our goal is then to learn a decision boundary in the feature space to differentiate between positive (with class label 1) and negative (with class label −1) training instances. This becomes a typical binary classification problem and any classification algorithm can be used for this task. We will use SVM to learn a decision boundary for this task in Section 3.2. For a given large TSP instance where we do not know the optimal route, the trained model can then be used to predict a class label for each edge in the graph. By removing the edges that are predicted to be −1, we have a reduced sparse graph, which is hopefully much easier for an existing solution algorithm to solve. The main steps of our MLPR method for TSP are summarized in Algorithm 1.
Algorithm 1 MLPR FOR TSP 1: Solve selected easy TSP instances to optimality; 2: Assign a class label 1 to the edges in the optimal route and −1 to others; 3: Extract features to characterize each edge (training instance); 4: Train a machine learning model to separate positive and negative edges; 5: Predict a class label for each edge on an unseen test graph (where its optimal solution is unknown) using the trained model, and remove negative edges.

Extracting Features to Characterize Each Edge
We extract four features directly computed from graph data and two statistical measures computed from stochastic samples of feasible solutions to characterize each edge (training instance).

Graph Features
As the objective of TSP is to search for a route with minimum costs, the edge cost is an important feature related to the objective value. Considering a TSP instance G(V, E,C) with n cities, we design four graph features to describe each edge e i, j , i, j = 1, · · · , n, based on the edge costs C: The first feature computes the difference between the edge cost of e i, j (c i, j ) and the minimum edge cost that originates from vertex i, while the second feature computes the difference between c i, j and the minimum edge cost that ends in vertex j. The third and fourth features are computed based on the mean edge costs connected to vertex i or j. We normalize the four features by the difference between the maximum and minimum edge costs that connect to vertex i or j. These graph features only capture local characteristics of an edge. In the next subsection, we will describe two statistical measures to capture certain global features for each edge.

Statistical Measures
The statistical measures aim to quantify the likelihood of each edge belonging to an optimal route based on randomly generated samples of feasible routes. As TSP has the backbone structure (Kilby et al., 2005), it is possible to identify the edges shared between an optimal route and high-quality routes.
Randomly generating a feasible route for TSP is very simple. Supposing the vertices (cities) are labelled from 1 to n, any random permutation (P) of integers from 1 to n is a feasible route for visiting each city. We generate m random feasible routes {P 1 , P 2 , · · · , P m }, and compute the corresponding objective values {y 1 , y 2 , · · · , y m }. The time complexity of sampling is Θ(mn), simply because a random permutation of n elements costs Θ(n).
To define our statistical measures, we introduce a binary string x x x k to represent the k th sample route P k , where x k i, j = 1 means the edge e i, j is in the k th sample; otherwise it is not. The first statistical measure is computed from the ranking of sample routes. We sort the sample routes based on their objective values in ascending order, and use r k to denote the ranking of the k th sample. The ranking-based measure for edge e i, j is defined as where i, j = 1, · · · , n. The edges that frequently appear in high-quality sample routes have a large rankingbased score, and are more likely to be part of an optimal route. We then normalize each ranking-based score by dividing the maximum ranking-based score in a graph This normalization avoids a large-valued feature dominating a classification task.
The second statistical measure we have developed is a correlation-based measure, that computes the Pearson correlation coefficient between each variable x i, j and objective values across the sample routes: wherex i, j = ∑ m k=1 x k i, j /m, andȳ = ∑ m k=1 y k /m. As TSP is a minimization problem, edges that are highly negatively correlated with the objective values are likely to be in an optimal route. Similarly, we normalize the correlation-based score by the minimum correlation value in a graph: Directly computing these two statistical measures from the binary string representation x costs Θ(mn 2 ) in both time and space complexity. In (Sun et al., 2019), we have introduced an efficient method based on set representation, i.e., permutation P here. Because x k i, j are binary variables, we can simplify the calculation of Pearson correlation coefficient using the following two equalities: where Algorithm 2 STATISTICAL MEASURES(P, Y , m, n) 1: Sort the samples in P based on objective value Y ; use r k to denote the ranking of k th sample P k ; 2: Compute mean objective value:ȳ ← ∑ m k=1 y k /m; 3: Compute objective difference: y d ← ∑ m k=1 (y k −ȳ); 4: Compute objective variance: σ y ← ∑ m k=1 (y k −ȳ) 2 ; 5: Initialize f r ,x i, j and s 1 i, j to 0, for each e i, j ∈ E; 6: for k from 1 to m do 7: for idx from 1 to n do 8: if idx < n then 9: ; 15: for i from 1 to n do 16: for j from 1 to n and j = i do 17: 20: return f r and f c .
The proof of these two equalities can be found in (Sun et al., 2019). Having this simplification, we can compute the two statistical measures in Θ(mn + n 2 ) for both time and space complexity using Algorithm 2. The main idea is to iterate through the edges in each sample route P to accumulate our ranking-based measure f r (e i, j ),x i, j and s 1 i, j , i.e., line 6 to 14 in Algorithm 2. Our correlation-based measure f c (e i, j ) can then be easily computed based onx i, j and s 1 i, j . In practice, the sample size m should be larger than n; otherwise there will be some edges that are never sampled. Considering a symmetric TSP instance with n cities, the number of edges in the undirected complete graph is n(n − 1)/2. The total number of edges in m sample routes is mn. Thus each edge is expected to be sampled 2m/(n − 1) times. In our experiments, we set m = 100n.
It is noteworthy that the two statistical measures described here can be directly used as a problem reduction technique, e.g., we can remove the edges that are positively correlated with objective values from a graph. However, as we have shown in (Sun et al., 2019), the machine learning approach which takes these two statistical measures as features outperforms a single statistical measure for problem reduction. Thus we will simply use both statistical measures as inputs to our machine learning model in this paper.

Support Vector Machine Classification
In our training set, the number of positive training instances is much smaller than that of negative instances. Considering a symmetric TSP instance with n cities, the number of edges in an optimal route is n, and the total number of edges is n(n − 1)/2. Thus the ratio between positive and negative edges is 2 : (n − 3). The standard SVM formulation tends to classify negative training instances better than the positive instances, because there are more negative training instances. However, misclassifying a positive instance is much harmful than misclassifying a negative instance. If a positive instance is misclassified, the reduced optimization problem no longer captures the original optimal solution. On the other hand, misclassifying a negative instance only results in a slight increase of the reduced problem size. In this sense, we will use cost-sensitive SVM (see Section 2.3) and penalize misclassifying positive instances more by using a larger regularization parameter r + in Eq. (14), in contrast to that of negative instances r − . In our experiments, we will set r − = 1 and r + = ε m n −1 /n 1 , where n −1 and n 1 are the number of negative and positive instances in our training set, and ε m controls the penalty for misclassifying positive instances. The term n −1 /n 1 balances the number of positive and negative instances in our training set.
We will consider two types of SVM in our experiments, linear SVM (solving primal optimization problem) and non-linear SVM with the RBF kernel (solving dual optimization problem); see Section 2.3 for more details. The classification accuracy of kernel SVM is usually higher than that of linear SVM. We will use the SMO-type (Sequential Minimal Optimization) decomposition method (Fan et al., 2005) implemented in the LIBSVM library (Chang & Lin, 2011) to solve the dual optimization problem of L1-SVM. However when the number of training instances is too large (e.g, millions of instances), solving the dual problem is computationally very slow. In this case, we will solve the primal optimization problem of linear L2-SVM using the trust region Newton method  implemented in the LIBLINEAR library (Fan et al., 2008).

Experiments
In this section, we use simulation experiments to investigate the robustness of our MLPR model. We will consider three scenarios where training and test instances are different, and explore the corresponding generalization errors. Specifically in Section 4.1, we train our MLPR model using one category of TSP instances and test it on another with different problem characteristics. In Section 4.2, we train MLPR using small randomly generated TSP instances, and test it on larger randomly generated TSP instances. In Section 4.3, we train our MLPR model on symmetric TSP instances, and test it on three TSP variants. In the last subsection, we investigate whether our MLPR method can be used as a preprocessing technique to boost the performance of a generic solver -CPLEX. Our source codes are implemented in C++, and complied with GCC/7.3.0-2.30. All our experiments are conducted on a high performance computing server with Intel(R) Xeon(R) Gold 6154 CPUs @ 3.00GHz and 21GB RAM.
Because the dimensionality of these instances is small (i.e., 100), we can quickly solve these instances to optimality using the Concorde solver (Applegate et al., 2006a). This enables us to systematically evaluate the generalization capability of our MLPR model to instances in different categories. We train a MLPR Table 1: The average optimality gap (%) generated by our MLPR model when trained on one category of TSP instances and tested on another. The 7 categories of instances are labelled as I 1 , · · · , I 7 . The MLPR model trained on category I j is denoted as MLPR-I j . The best optimality gap generated for each test category is in bold.

Models
I 1 I 2 I 3 I 4 I 5 I 6 I 7  model using the first 50 instances in one category, and test the trained model on the remaining unseen 140 instances in that category as well as the instances from other 6 categories. For each test instance, we apply the trained model to reduce the problem size (i.e., pruning some edges in the complete graph), and solve the reduced problem to optimality using Concorde. We compute an optimality gap by comparing the optimal solutions generated in the reduced and original problems. As our MLPR model uses statistical features computed from random samples, the reduced problem generated for a test instance might be slightly different if we use a different random seed. Therefore, we repeat the random sampling and problem reduction process 25 times to alleviate randomness.

MLPR
As the size of each training set is not too large, we use L1-SVM with RBF kernel to train our MLPR model. The computational time required for training one model is less than 30 minutes, and the predicting time for one test instance is around 5 seconds. The kernel parameter γ k is set to the default value used in the LIBSVM library: γ k = 1/n f , where n f is the number of features. We have tested multiple values for the penalty parameter ε m in (Sun et al., 2019), and found that ε m = 10 works reasonably well across a wide range of problem instances. Thus we simply set ε m = 10 in this subsection.

Results
The average optimality gap generated by our MLPR model and remaining problem size after reduction are presented in Table 1 and 2 respectively. These two tables should be read in conjunction as, in an ideal world, we would want to see both small reduced problems and small optimality gaps. What the results show in practice is that, unsurprisingly, those MLPR models that retain a larger fraction of the original problem also tend to produce smaller gaps. Hence the best optimality gaps in Table 1 do not simply occur on the diagonal as one might first expect. As a compromise between these two aims of minimising gaps and size, we observe that using randomly generated TSP instances (I 7 ) as a training set, our MLPR model performs reasonably well across all 7 categories of TSP instances. The MLPR-I 7 model prunes about 85% of edges, and overall achieves 0.44% optimality gap. Further, the TSP instances that are hard for CLK or LKCC to solve (e.g., I 2 and I 6 ) are also hard for our MLPR model. The optimality gaps generated by our MLPR models for hard TSP instances are larger than those for easy TSP instances. Third, MLPR trained on easy TSP instances (e.g., MLPR-I 1 or MLPR-I 5 ) prunes too many edges for hard instances (e.g., I 2 or I 6 ), resulting in a large optimality gap. On the other hand, MLPR trained on hard TSP instances prunes too few edges for easy instances, resulting in a large remaining problem size. For example, both MLPR-I 2 and MLPR-I 3 achieve 0 optimality gap for easy instances in I 1 , but MLPR-I 2 prunes 8.5% few edges than MLPR-I 3 . Last, we observe that when training and test instances are from the same category, the generated optimality gap is less than 0.2%, and the amount of pruned edges is reasonable -in fact the results are always non-dominated in the Pareto sense, with all other models producing worse gaps or larger remaining problems. We use the tool from the MATILDA library to visualize the performance of our MLPR models. Here each TSP instance (instead of an edge) is mapped to a point in a 6-dimensional feature space. The feature values are taken from the MATILDA library and are normalized to the range of 0 to 1. We then apply principal component analysis to reduce the 6-dimensional feature space to 2-dimensional by selecting the first two principal components (Z 1 and Z 2 ). The distribution of the TSP instances in the 2-dimensional space, spanned by Z 1 and Z 2 , is shown in Figure 1a. For better visualization, we only plot the last 50 TSP instances from each category, and the instances from different categories are in different colors. We can see that instances from different categories are well separated in the feature space. Note that we only visualize five categories of TSP instances, because the feature data for the other two categories is not available in MATILDA.
We also plot the optimality gap (%) generated by our MLPR models for these instances when trained on one category of instances at a time in Figure 1b-1f. The optimality gap is indicated by the color of the dots (blue is small and red is large). We can observe that the MLPR model trained on easy TSP instances (MLPR-I 1 or MLPR-I 5 ) does not perform well on hard instances (I 2 or I 6 ). There appears to be a strong correlation between the size of the gap and the distance from the training data for these instances. The MLPR model trained on hard instances is able to generate a small optimality gap for any instance considered, though this comes at a price of an increased problem size after the reduction.

Setting
In this subsection, we explore the generalization error of our MLPR model in terms of problem size. We train our MLPR model using 190 randomly generated TSP instances (I 7 ) from the MATILDA library with dimension 100, and test it on larger randomly generated TSP instances. Each TSP instance in I 7 of MATILDA is created by randomly generating 100 pairs of integer coordinates between 0 and 400. We use the same method to create test instances with different number of cities, i.e., 200, 500, 800, 1100, 1400, 1700, and 2000. This results in 7 categories, each with 190 randomly generated TSP instances. Because each edge (instead of a complete graph) is a training instance, our training set size is close to 1 million. Thus we use linear SVM with ε m = 10 to train our MLPR model, by solving the primal problem to gain computational efficiency. We apply the trained model to reduce problem size for each test TSP instance, and solve the reduced problem as well as the original problem to optimality using Concorde.

Results
Surprisingly, the optimality gap generated by our MLPR model for each test category is always zero. As the dimension of the test TSP instances increases, our MLPR model tends to prune slightly more edges, as shown in Figure 2. This suggests that our MLPR model trained on small randomly generated TSP instances generalize well to larger randomly generated TSP instances.

Setting
In this subsection, we train our MLPR model on the symmetric TSP and test it on other TSP variants. We take the first 50 LKCC-hard instances (I 6 ) from MATILDA as our training set, because it overall generates the smallest optimality gap for symmetric TSP instances with different problem characteristics according to Table 1. We train our MLPR model using SVM with an RBF kernel and test two penalty parameter values ε m = 10 and ε m = 100. We gradually differentiate the test instances away from the training instances, by using 1) symmetric TSP, 2) asymmetric TSP, and 3) sequential ordering problem (SOP) instances in testing. The asymmetric TSP is a TSP variant that allows the distance matrix to be asymmetric; and SOP is a variant that further considers precedence constraints in the order of visiting cities. These test instances are all from the TSP library (Reinelt, 1991). For symmetric TSP, we use 19 instances for which the number of cities is in between 100 and 200. For asymmetric TSP and SOP, we use the easy instances that can be solved to optimality by CPLEX with 8 CPUs in 1000 seconds. We use the trained MLPR model to reduce problem size for each test instance, and solve the original and reduced problems to optimality by exact solvers (Concorde for symmetric TSP, and CPLEX for asymmetric TSP and SOP instances). The mixed integer formulation and random sampling method used for SOP are presented in Appendix A and B. As before, the random sampling and problem reduction process is repeated 25 times to allow for randomness.

Results
The optimality gap generated by our MLPR model and remaining problem size after reduction for each test instance are presented in Figure 3, and the average statistics across each problem type are presented in Figure 4. These results show that the MLPR model trained on symmetric TSP instances makes useful predictions about which variables can be eliminated without significantly impacting solution quality when testing on instances from asymmetric TSP and SOP. We also observe that as we gradually move the test instances away from the training instances, our MLPR model becomes less accurate, resulting in a larger optimality gap. On the other hand, it also becomes less confident at pruning edges, resulting in a larger remaining problem size. When we use a larger penalty ε m = 100, our MLPR model prunes less edges so the optimality gap generated is smaller.
When the training and test instances are both from symmetric TSP, our MLPR model performs well; it prunes on average 75.89% of edges but still captures a near-optimal solution in the reduced problem, which is within 0.56% from the optimal solution.
When training on symmetric TSP and testing on asymmetric TSP instances, the optimality gaps generated by our MLPR model (with ε m = 100) are all less than 1%, except for two instances rbg403 and rbg443. We then use two closely-related asymmetric TSP instances rbg323 and rbg358 as training set, and investigate whether the performance of our MLPR model can be improved for rbg403 and rbg443. The results show that by using asymmetric TSP as training instances, the optimality gap generated by our MLPR model (with ε m = 10) for rbg403 and rbg443 can be significantly reduced from 10.39% to 1.2% and from 11.91 % to 1.25% respectively.
When training on symmetric TSP and testing on SOP instances, our MLPR model does not perform very well on some instances e.g., ESC12, br17.10 and rbg109a. However we show that by using the other 6 SOP instances as training set, the optimality gaps generated by our MLPR model (with ε m = 10) for ESC12, br17.10 and rbg109a can be significantly reduced to 0.65%, 0.00% and 0.23% respectively.

Setting
In this subsection, we investigate whether our MLPR model can be used as a preprocessing technique to boost the performance of a generic solver -CPLEX. Because the symmetric TSP can be efficiently solved by a specialized solver -Concorde, and an asymmetric TSP can be easily converted to a symmetric TSP (Jonker & Volgenant, 1983), we mainly focus on solving SOPs here.
We take the 9 easy SOP instances used in Section 4.3 as our training set, and train a machine learning model using SVM with RBF kernel (ε m = 10), which costs a few minutes. We then apply the trained model to reduce problem size for hard SOP instances from the TSP library. We select 24 hard SOP instances for testing, which have not been proved to optimality according to the bounds published in the TSP library. The reduction time for the largest test instance (with 378 cites) is around 25 seconds. We use CPLEX with 8 CPUs to solve the reduced problem, compared to directly solving the original problem. The cutoff time is set to 1000 seconds for both, and the reduction time used by our MLPR model is counted as part of the cutoff time. Note that we also parallelize the process of generating random samples and computing statistical features for our MLPR model. The MIP emphasis parameter is set to "HIDDENFEAS" and MIP search method is set to "TRADITIONAL" for CPLEX. The experiments are again repeated 25 times with different random seeds. We report the successful rate of finding any feasible solution for the 25 independent runs and compute an average optimality gap only for the successful runs.

Results
The optimization results of CPLEX and CPLEX-MLPR when used to solve the SOP instances are presented in Table 3. We observe that by using our MLPR method as a preprocessing technique, CPLEX can generally find a comparable or better solution within the cutoff time, especially for hard instances. For 13 out of 24 instances, CPLEX fails to find a feasible solution in any of the 25 runs, i.e., the successful rate is 0. In contrast, CPLEX-MLPR can always find a feasible solution. This is partially because we feed the best feasible solution from sampling to CPLEX as a warm start.
It is important to note that when using our MLPR method as a preprocessing technique, the generic solver CPLEX can find a better primal solution that improves the best objective value published in the TSP library for 4 instances (highlighted in bold in Table 3). The best objective values found by CPLEX-MLPR for some instances are much larger than the lower bounds published in the TSP library, resulting in a large optimality gap (gap b and gap m ). However this is mainly due to that CPLEX is ineffective at solving these instances, which can be inferred from the fact that the lower bounds produced by CPLEX-MLPR are usually much smaller than the lower bounds published in the TSP library. Thus CPLEX-MLPR can potentially find a better primal solution if given enough computational time. We have also verified that the lower bounds produced by CPLEX-MLPR are less than the best objective values published in the TSP library for all of these instances. This indicates our problem reduction method MLPR does not introduce a noticeable error when pruning edges. Table 3: The optimization results (in percentage) of CPLEX and CPLEX-MLPR when used to solve the hard SOP instances within the cutoff time 1000 seconds. The columns gap b and gap m compute the best and average optimality gaps of the feasible solutions found by CPLEX or CPLEX-MLPR comparing to the lower bounds published in the TSP library; while gap s computes an optimality gap based on the lower bounds generated by CPLEX-MLPR; r is the successful rate of finding a feasible solution; p is the percentage of remaining problem size after reduction. The last column presents the optimality gap published in the TSP library. We highlight the 4 instances for which the best optimality gap generated by CPLEX-MLPR is better than the gap given in the TSP library.

Conclusion
In this research, we have applied machine learning techniques to reduce problem size for combinatorial optimization, which can be used as an effective preprocessing step to improve the performance of existing solution algorithms. We have adapted the machine learning model that we originally developed in (Sun et al., 2019) to prune the search space for the travelling salesman problem (TSP). We empirically investigated the generalization error of our machine learning model when training and test (unseen) instances have different instance characteristics, sizes or are from different TSP variants. Our experimental results showed that our model generalized reasonably well to a wide range of instances with different characteristics or sizes. In general, when training and test instances are from the same TSP variant, the generalization error of our model is small; and this generalization error naturally increases when testing our model on TSP variants that become increasingly dissimilar to the training instances. Solving a completely different problem variant than in the training set is a fairly extreme case of mismatch between training and application of the model. Hence this shows that our approach is likely to be quite robust to the more typical changes seen in practice when applying a model to, e.g., a business application where the problem characteristics slowly drift over time. In future, we plan to develop a more generic model that does not require re-training when applied to solve a class of problems. Another possible direction for future work would be to apply our problem reduction method to real-world problems where a fast computation of a high-quality solution is desired, such as trip planning.