Year: 2021

Abstract:
Deep reinforcement learning (RL) has proved to be a competitive heuristic for solving small-sized instances of traveling salesman problems (TSP), but its performance on larger-sized instances is insufficient. Since training on large instances is impractical, we design a novel deep RL approach with a focus on generalizability. Our proposition consisting of a simple deep learning architecture that learns with novel RL training techniques, exploits two main ideas. First, we exploit equivariance to facilitate training. Second, we interleave efficient local search heuristics with the usual RL training to smooth the value landscape. In order to validate the whole approach, we empirically evaluate our proposition on random and realistic TSP problems against relevant state-of-the-art deep RL methods. Moreover, we present an ablation study to understand the contribution of each of its components.

heading: 1 INTRODUCTION
text: The traveling salesman problem (TSP) has numerous applications from supply chain management (Snyder & Shen, 2019) to electronic design automation (Gerez, 1999) or bioinformatics (Jones & Pevzner, 2004). As an NP-hard problem, solving large-sized problem instances is generally intractable. Deep reinforcement learning (RL)-based heuristic solvers have been demonstrated (Bello et al., 2016; Dai et al., 2017; Kool et al., 2019; Ma et al., 2019) as being able to provide competitive solutions while being fast, which is crucial in many domains such as logistics for real-time operations. However, this ability has only been demonstrated on small-sized problems where RLbased solvers are usually trained and tested on small instances. When evaluated instead on larger instances, such solvers perform poorly. Since the computational cost of training on large instances is prohibitive, increasing the size of the training instances is not a practical option for obtaining efficient solvers for larger instances. To overcome this limitation, this paper investigates techniques to increase the generalization capability of deep RL solvers, which would allow to train on small instances and then solve larger ones.
Generalization in deep RL (François-Lavet et al., 2018) can be improved by adjusting the following three components: input/representation space, objective function, and learning algorithm. We propose to achieve this by exploiting equivariance, local search, and stochastic curriculum learning, which we explain next. Although our approach is demonstrated on TSP problems, we believe it can be naturally extended to other combinatorial optimization problems, but also to more classic RL domains. We leave this for future work.
Problems with spatial information, such as TSP problems, enjoy many spatial symmetries which can be exploited for RL training and generalization. An RL solver is invariant with respect to some symmetry, if its outputs remain unchanged for symmetric inputs. For instance, translating all the city positions leave optimal solutions intact. More generally, an RL solver is equivariant, if its outputs are also transformed with a corresponding symmetry for symmetric inputs. For instance, if cities are permuted, a corresponding permutation is required on the outputs of the RL solver. Invariance and equivariance of the solver with respect to some symmetries allows a smaller input space to be considered during training and more abstract representation insensitive to those symmetries to be learned while the trained solver still covers the whole input space.
Previous work has considered using local search as a simple additional step to improve the feasible solution output by a RL-based solver. In contrast to most previous work, we further employ local
search as a tool to smooth the objective function optimized during RL training. To that aim, we interleave RL training with local search using the improved solution provided by the latter to train the RL solver via a modified policy gradient. Moreover, we propose a novel simple baseline called the policy rollout baseline to reduce the variance of its estimation.
Curriculum learning (Soviany et al., 2021) can accelerate training, but also improve generalization (Weinshall et al., 2018). By stochastically adjusting the difficulty of the training instances, as described by instance sizes, the RL-based solver can learn more easily and faster. Since the solver sees various sizes of instances, this may prevent it to overfit to one particular instance size.
The architecture of our model includes a graph neural network (GNN), a multi-layer perceptron (MLP), and an attention mechanism. Due to all the used techniques, we name our model as eMAGIC (e for equivariance, M for MLP, A for Attention , G for Graph neural network, I for Interleaved local search, and C for Stochastic Curriculum learning). We demonstrate that it can be trained on small random instances (up to 50 cities) and perform competitively on larger random or realistic instances (up to 1000 cities). The only learning methods that can perform better than ours require learning on the instance to be solved (see Section 2 for a discussion).

heading: 2 RELATED WORK
text: Research on exploiting deep learning (Vinyals et al., 2015; Li et al., 2018; Joshi et al., 2019a; Prates et al., 2019; Xing & Tu, 2020) or RL (Bello et al., 2016; Dai et al., 2017; Kool et al., 2019; Deudon et al., 2018; Ma et al., 2019; Da Costa et al., 2020; Zheng et al., 2021; Kwon et al., 2021; Wu et al., 2021b) to design heuristics for solving TSP problems has become very active. The current achievements demonstrate the potential of machine-learning-based approaches, but also reveal their usual limitation on solving larger-size instances, which has prompted much research on generalization (Joshi et al., 2019b; Fu et al., 2021; Wu et al., 2021b). Some experimental work (Joshi et al., 2019b) suggests that deep RL may provide more generalizable solvers than supervised learning, which is a further motivation for our own work. Another advantage of deep RL is that it does not require optimal solvers to train. We discuss next the related work that exploits equivariance or local search like our method. To the best of our knowledge, no other work uses stochastic curriculum learning for designing an RL-based solver for TSP, although Lisicki et al. (2020) evaluated deterministic curriculum learning strategies on small TSP instances. For space reasons, our discussion below emphasizes the approaches using deep RL as a constructive heuristic (which builds a solution iteratively).
Invariance and equivariance have been important properties to exploit for designing powerful deep learning architectures (Bengio & Lecun, 1997; Gens & Domingos, 2014; Cohen & Welling, 2016). They have also inspired some of the previous work on solving TSP problems, although they were often not explicitly discussed like in our work. For instance, Deudon et al. (2018) use principal component analysis to exploit rotation invariance as a single preprocessing step. In contrast to their work, we compose various preprocessing transformations to be applied at every solving step. To the best of our knowledge, we are the first to propose such technique. In addition, permutation invariance is the motivation for using attention models (Kool et al., 2019; Deudon et al., 2018) or GNNs (Ma et al., 2019) in RL solvers. Contrary to Deudon et al. (2018), Kool et al. (2019) select the next city to visit based on the first and last visited cities while Ma et al. (2019) propose to use relative positions for translation invariance. Our proposition combines those two ideas, but compared to Kool et al. (2019), our model is based on GNN, while compared to Ma et al. (2019), it has a simpler architecture that does not require an LSTM. For vehicle routing problems, Peng et al. (2020) propose to remove visited cities from the RL agent’s input, a technique we also apply to TSP. In contrast to their work, our GNN model is simpler and we investigate generalization from small instances to very large instances. To summarize, compared to all previous work, our paper investigates in a more systematic fashion the exploitation of equivariance to help training and improve generalization.
Several previous studies combine deep learning or RL with various local search techniques to find better solutions: 2-opt (Deudon et al., 2018; Ma et al., 2019; Da Costa et al., 2020; Wu et al., 2021b), k-opt (Zheng et al., 2021), or Monte-Carlo tree search (Fu et al., 2021). In contrast to these methods, we use a combined local search technique, which applies several heuristics in an efficient way. Moreover, to the best of our knowledge, no previous work considers interleaving local search with policy gradient updates to obtain a more synergetic final method for solving TSP. For the bin packing problem, Cai et al. (2019) use RL to learn a perturbative heuristics to provide a good initial
solution to a heuristic optimizer. Our interleaved training process is based on a similar idea, but our RL agent learns a constructive heuristics and we provide an intuitive motivation why this approach may work via our smoothed policy gradient.
All these machine-learning-based methods can be categorized as learning from a batch of instances or as directly learning on the instance to be solved. Like our work, most methods are part of the first category. However, some recent work belongs to the second category (Zheng et al., 2021) or are hybrid (Fu et al., 2021). For instance, Zheng et al. (2021) apply RL to make the LKH algorithm (Helsgaun, 2017), a classic efficient heuristic for TSP, adaptive to the instance to be solved. In Fu et al. (2021)’s method, evaluations learned in a supervised way from a batch of instances are used in Monte-Carlo tree search so that solutions can adaptively be found for the current instance. Being adaptive to the current instance can undoubtedly boost the performance of the solver. We leave for future work the improvement of our solver to become more adaptive to the instance to be solved.

heading: 3 PRELIMINARIES
text: Following previous work, we focus on the symmetric 2D Euclidean TSP, which we recall below. We then explain how a TSP instance can be tackled with RL. We first provide some notations: For any positive integerN ∈ N, [N ] denotes the set {1, 2, . . . , N}. Vectors and matrices are denoted in bold (e.g., x or X). For a set of subscripts I , XI denotes the matrix formed by the rows of X whose indices are in I .
Traveling Salesman Problem A symmetric 2D Euclidean TSP instance is described by a set of N cities (identified to the set [N ]) and their coordinates in R2. The goal in this problem is to find the shortest tour that visits each city exactly once. The coordinates of city i are denoted xi ∈ R2. The matrix whose rows correspond to city coordinates is denoted X ∈ RN×2. A feasible solution (i.e., tour) of a TSP instance is a permutation σ over [N ] with length equal to:
Lσ(X) = N∑ t=1 ‖xσ(t) − xσ(t+1)||2 (1)
where ‖·‖2 denotes the `2-norm, σ(t) ∈ [N ] is the t-th visited city in tour σ, and by abuse of notation, σ(N + 1) denotes σ(1). Therefore, solving a TSP instance consists in finding the permutation σ that minimizes the tour length Lσ(X) defined in Equation (1). Since TSP solutions are invariant to scaling, we assume that the city coordinates are in the square [0, 1]2. In Appendix A, we recall some heuristics (e.g., insertion, k-opt) that have been proposed for solving TSP.
RL as a Constructive Heuristic RL can be used to construct a TSP tour σ sequentially. Intuitively, at iteration t ∈ [N ], an RL solver (i.e., policy) selects the next unvisited city σ(t) to visit based on the current partial tour and the description of the TSP instance (i.e., coordinates of cities). Therefore, this RL problem corresponds to a repeated N -horizon sequential decision-making problem. As noticed by Kool et al. (2019), the decision for the next city to visit only depends on the description of the TSP instance and the first and last visited cities. In addition, we update the description of the TSP instance by removing the coordinates of visited cities. Surprisingly, to the best of our knowledge, no previous work exploits this simplified RL model.
Formally, in the RL language, at time step t ∈ [N ], an action at represents the next city to visit, i.e., at = σ(t). Let I1 = [N ] and It+1 = [N ]\{σ(1), . . . , σ(t)} be the set of remaining cities after t cities have been already visited. Moreover, let J1 = [N ] and Jt = It ∪ {σ(1), σ(t)} be the set of unvisited cities in addition of the first and last visited cities (σ(1) and σ(t)). Therefore, at ∈ It for all t ∈ [N ]. A state st can be represented by a matrix XJt with flags indicating the first and last visited cities. This matrix includes the coordinates of unvisited cities (XIt ) in addition to the coordinates of the first and last visited cities (xσ(1) and xσ(t)). Note that at t = 1, no city has been chosen yet, so the initial state s1 only contains the list of city coordinates. A state sN contains the coordinates of the unvisited cities and those of the first and last visited cities. After choosing the last city to visit in state sN , the tour σ is completely generated. The immediate reward r(st, at) for an action at in a state st can be defined as the negative length between the last visited city and the next
chosen city, since we want the tour length to be small:
r (st, at) =  0 for t = 1 −||xσ(t) − xσ(t−1)||2 for t = 2, . . . , N − 1 −||xσ(N) − xσ(N−1)||2 − ||xσ(1) − xσ(N)||2 for t = N
(2)
After choosing the first city a1 = σ(1), the reward is zero since no length can be computed. After the final action aN , an additional reward −||xσ(1) − xσ(N)||2 is added to complete the tour length. In deep RL, a policy πθ is represented as a neural network parametrized by θ. The goal is then to find θ∗ that maximizes the objective function J(θ):
θ? = argmax θ J(θ) = argmax θ Eτ∼pθ(τ) [ N∑ t=1 rt ] , (3)
where for all t ∈ [N ], rt = r(st, at), τ = (s1, a1, s2, a2, . . . , sN , aN ) is a complete trajectory, and pθ is the probability over trajectories induced by policy πθ. This objective function J(θ) can be optimized by a policy gradient (Williams, 1992) or actor-critic method (Sutton & Barto, 1998).

heading: 4 EQUIVARIANT MODEL
text: In this section, we present several techniques to exploit invariances and equivariances in the design of an RL solver. Formally, a mapping f : A → B from a set A to a set B is invariant with respect to a symmetry ρA : A→ A iff f(x) = f(ρA(x)) for any x ∈ A. More generally, given a symmetry that acts on A with ρA : A→ A and on B with ρB : B → B, a mapping f : A→ B is equivariant with respect to this symmetry iff ρB(f(x)) = f(ρA(x)) for any x ∈ A. This definition shows that invariance is a special case of equivariance when ρB is the identity function. Intuitively, equivariance for an RL solver (resp. its value function) means that if a transformation is applied to its input, its output (resp. its value) can be recovered by a corresponding transformation.
In the remaining of the paper, for simplicity, we often use equivariance to refer to both equivariance and invariance, since the latter is a special case of the former. Equivariance can be exploited in RL in various ways. We consider some equivariant preprocessing methods on the description of TSP instances to standardize the kinds of instances the solver is trained on and evaluated on. In addition, we propose a simple deep learning model for which we apply some other equivariant preprocessing methods on its inputs to further reduce the input space.

heading: 4.1 EQUIVARIANT PREPROCESSING OF TSP INSTANCE
text: An RL solver should be invariant with respect to any Euclidean symmetry (rotation, reflection, translations and their compositions) and to any positive scaling transformation applied on city positions. To enforce these invariances, we can apply these transformations to preprocess the inputs of the solver such that the transformed inputs are always in a standard form. Doing so allows the solver to be trained on more similar inputs.
Concretely, for rotation invariance, we rotate and scale the city positions such that they are mostly distributed along the first diagonal of the square [0, 1]2. This can be achieved by performing a principal component analysis, rotating the first found axis by 45° anti-clockwise and scaling to fit the cities in the [0, 1]2 square. This transformation, which is slightly different from (Deudon et al., 2018), allows the cities to be as spread as possible in [0, 1]2. For scaling and translation invariance, we apply a scaling and translation transformation to the city positions such that there are a maximum number of cities (i.e., 2 or 3 depending on configuration) on the border of the square [0, 1]2. For reflections, we only consider horizontal, vertical, and diagonal flips (with respect to the coordinate axes) for simplicity. Cities are reflected such that a majority of them are in a fixed chosen region.
Since the symmetries can be composed, these preprocessing methods can be applied sequentially. Theoretically, it would be beneficial to apply all of them in combination, however this has a computational cost. It is therefore more effective to only use a selection of the most effective ones. Moreover, these preprocessing methods can be applied either once on the initial TSP instance (X), or at each solving step t on the remaining cities (XIt ). Theoretically, performing these preprocessings iteratively should help most, since the RL solver is only trained on standardized inputs.
This is confirmed by our experimental analysis. We find out that the combination that provides the best results are rotation followed by scaling and translation, which will be used in the experimental evaluation of our method. However, our overall approach is generic and could include any other symmetries for which equivariance would hold for the RL solver. In Appendix F, we present the details of the evaluation of the different preprocessing methods.
In the remaining, we assume that the set of equivariant transformations is fixed. For any set of indices I ⊆ [N ], we denote the matrix of positionsXI after preprocessing by X̃I and any coordinates x after preprocessing by x̃. Note that the preprocessing step (e.g., if it includes the rotation transformation described above) may depend on the initial matrix XI , but we prefer not to reflect it in the notations to keep them simple.

heading: 4.2 EQUIVARIANT MODEL
text: Our proposed model, which represents the RL solver, has an encoder-decoder architecture (see Figure 1). In our model, the first city to visit is fixed arbitrarily, since the construction of the tour should be invariant to the starting city. Therefore, the decisions of the solver starts at time step t ≥ 2. Our model is designed to take into account other equivariances that are known to hold in the problem. Invariance with respect to translation can be further exploited by considering relative positions with respect to the last visited city instead of the original absolute positions, as suggested by Ma et al. (2019). Therefore, we define the input of our RL solver to only include two pieces of information. First, the information about the current partial tour now only needs the relative preprocessed position of the first visited city (~xσ(1) = x̃σ(1) − x̃σ(t)). Since the last visited city is always represented by the origin, it can be dropped. Second, the information about the remaining TSP instance corresponds to the relative preprocessed positions of the unvisited cities (~xσ(n) = x̃σ(n) − x̃σ(t−1) for n ∈ It) and of the first and last visited cities (~xσ(1) and the origin ~xσ(t) = 0). We denote the matrix of those relative preprocessed positions by −→ XJt , i.e., the matrix whose rows are the rows of X̃Jt (obtained from XJt ) minus x̃σ(t). The last visited city needs to be kept in this matrix, since it represents the labels of the nodes of the remaining graph on which the tour should be completed.
Encoder With H denoting the embedding dimension, the encoder of our model is composed of a graph neural network (GNN) (Battaglia et al., 2018), which computes an embedding X̂Jt ∈ R|Jt|×H of the relative city positions −→ XJt ∈ R|Jt|×2, and a multilayer perceptron (MLP), which computes an embedding x̂σ(1) ∈ RH of the relative position of the first visited city −→x σ(1) ∈ R2. The GNN encodes the information about the graph describing the remaining TSP problem. As a GNN, it is equivariant with respect to the order of its inputs and its outputs depend on the graph structure and information (i.e., city positions). The MLP encodes the information about the first city. Since we exploit the independence with respect to the visited cities between the first and last cities (not included), our model does not require any complex architecture, like LSTM (Hochreiter & Schmidhuber, 1997) or GNN, for encoding the positions of the visited cities, in contrast to most previous work. A simple MLP suffices since only the relative position of the first city is required.
Formally, the GNN computes its output X̂Jt ∈ R|Jt|×H from inputs −→ XJt ∈ R|Jt|×2 through nGNN layers, which allows each city to get information from its neighbors up to nGNN edges away:
X(0) = −→ XΘ(0) (4)
X(`) = λ ·X(`−1)Θ(`) + (1− λ) · F(`) ( X(`−1)
|Jt| − 1
) (5)
where Θ(0) ∈ R2×H and Θ(`) ∈ RH×H are learnable weights, X(`−1) ∈ R|Jt|×H is the input of the `th layer of the GNN for ` ∈ [nGNN], X(nGNN) = X̂Jt , F(`) : R|Jt|×H → R|Jt|×H is the aggregation function, which is implemented as a neural network, and λ ∈ [0, 1] is another trainable parameter. Visual illustration of the GNN is deferred to Appendix H for space reasons.
Decoder Once the embeddings X̂Jt and x̂σ(0) are computed, the probability of selecting a next city to visit is obtained via a standard attention mechanism (Bello et al., 2016) using X̂Jt and x̂σ(0) as the keys and query respectively. Formally, the decoder outputs a vector u ∈ RN expressed as:
uj = { −∞ ∀j ∈ It w · tanh ( X̂Jt,jΘg + Θmx̂ ) otherwise
(6)
where uj is the jth entry of the vector u, X̂Jt,j is the j th row of the matrix X̂Jt , Θg and Θm are trainable matrices with shape H ×H , w ∈ RN is another trainable weight vector. Then, a softmax transformation turns u into a probability distribution p = ( pj ) j∈[N ] over unvisited cities:
p = softmax(u) = ( euj∑N j=1 e uj ) j∈[N ]
(7)
where pj is the jth entry of the probability distribution p and uj is the jth entry of the vector u. Note that the probability of any visited city j ∈ It is zero since uj = −∞ in Equation (6).

heading: 5 ALGORITHM & TRAINING
text: We introduce three innovative training techniques that we apply during our training process: combined local search, smoothed policy gradient, and stochastic curriculum learning. Overall, these techniques can help train a fast and generalizable policy, which is able to generate tours that can easily be improved by local search. The overall algorithm1 can be found in appendix C.
Combined Local Search In contrast to previous work, which generally only considers 2-opt, we propose to use a combination of several (possibly random) local search methods (i.e., local insertion heuristic, random 2-opt, search 2-opt and search random 3-opt) to improve a tour generated by the RL solver. For space reasons, we defer their presentations in Appendix A. A local search method can heuristically improve a tour, but may get stuck in a local optimum. Using a combination of them alleviates this issue, since different heuristics usually have different local optima. Another important distinction with previous work is that we also use local search during training, not only testing.
Smoothed Policy Gradient For simplicity, we train our model with the REINFORCE algorithm (Williams, 1992), which leverages the policy gradient for optimization. However, instead of using the standard policy gradient, which is based on the value of the tour generated by the policy, we compute the policy gradient with the value of the tour improved by our combined local search. While standard RL training may yield policies whose outputs may not easily be improved by local search, our new definition directly trains the RL solver to find solutions that can be improved by local search, which allows RL and local search to have synergetic effects.
Intuitively, this new policy gradient amounts to smoothing the objective function J(θ) that is optimized. Recall J(θ) = −E[Lσ(X)], where the expectation is taken with respect to σ, which is a
1Our implementation takes advantage of GPU acceleration when possible. The source code will be shared after publication.
random variable corresponding to the tour generated by policy πθ. In our approach, this usual objective function is replaced by J+(θ) = −E[Lσ+(X)], where σ+ is a random variable corresponding to the improved tour obtained by our combined local search from σ. Therefore, this last expectation is taken with respect to the probability distributions generated by policy πθ and our combined local search. This objective function can be understood as J+(θ) = −E[minσ′∈N (σ) Lσ′(X)] where N (σ) is a neighborhood of σ defined by local search. This stochastic min operation has a smoothing effect on J(θ). That is why, we call the gradient of J+(θ) a smoothed policy gradient.
Formally, this novel policy gradient can be estimated on a batch of TSP instancesX(1), . . . ,X(B):
∇θJ+(θ) ≈ − 1
|B| |B|∑ b=1 ( N∑ t=2 ∇θ log πθ(a(b)t |s (b) t ) )( L σ (b) + (X(b))− l(b) ) , (8)
where σ(b) is the tour generated by the current policy πθ on instance X(b), σ (b) + is the improved tour obtained by our combined local search starting from σ(b), and l(b) = −Lσ(b)(X(b)) is a novel baseline, which we call policy rollout baseline, used to reduce the variance of the policy gradient estimation. It enjoys the nice property that it does not require additional calculations, since Lσ(b)(X
(b)) is computed when σ(b) is generated. In our experiments, the policy rollout baseline easily outperforms the previous greedy baselines (Ma et al., 2019; Kool et al., 2019).
Note that by construction, L σ (b) + (X(b)) ≤ Lσ(b)(X(b)) = −
∑N t=1 r (b) t . Therefore, the smoothed
policy gradient updates more if our combined local search can make more improvements upon this certain policy. For completeness, we provide more details in Appendix B.
Stochastic Curriculum Learning Curriculum Learning (CL) is widely used in machine learning (Soviany et al., 2021). Its basic principle is to control the increase of the difficulty of training instances. CL can speed up learning and improve generalization (Weinshall et al., 2018). We apply stochastic CL to train our model. Instead of a deterministic process, stochastic CL increases the difficulty according to a probability distribution. We choose the TSP size (i.e., number of cities) as a measure of difficulty for a TSP instance. We assume it to be in the range ofR = {10, 11, . . . , 50}. For each epoch e, we define the vector g(e) ∈ R41 (since |R| = 41) to be:
g (e) k = 1√ 2πσN exp − 12
( (k+10)−e
σN )2 , (9)
where g(e)k is the k-th entry of g (e) and hyperparameter σN is the standard deviation of this Gaussian density function. Then, g(e) is turned into a categorical distribution p(e) ∈ [0, 1]41 via a softmax:
p(e) = softmax(g(e)) (10)
The k-th entry of p(e) gives the probability of choosing a TSP instance of size (k + 10) at epoch e.

heading: 6 EXPERIMENTAL RESULTS
text: We present three sets of experiments. First, to validate the effectiveness of eMAGIC, we train our model on randomly-generated TSP instances (using stochastic CL with sizes up to 50), and test the model on other randomly-generated TSP instances (TSPn where size n = 20 up to 1000). Second, to further prove of its generalization capability, we directly evaluate models trained on random instances on realistic symmetric 2D Euclidean TSP instances with sizes range from 51 to 1002 in TSPLIB (Reinelt, 1991). Third, we conduct an ablation study to show the significance of every component of eMAGIC (i.e., equivariance, stochastic CL, policy rollout baseline, combined local search, and RL). We evaluate four versions of our model: eMAGIC(G), eMAGIC(S), eMAGIC(s1) and eMAGIC(s10) where G means the tour is generated greedily from the RL policy, while the other ones are based on random sampling and differ with respect to the number of times eMAGIC is applied (once for s1, 10 times for s10 and 100 for S). The details about the experimental settings and the used hyperparameters, which are the same for all experiments, are provided in Appendix D.
Performance on Randomly Generated TSP We compare the performance of our model with 12 other methods in Tables 1 and 2, which covers various types of TSP solvers including exact solvers,
heuristics, and learning-based approaches. Column 1 and 2 of both Tables 1 and 2 correspond to the method name and the method type, respectively. Columns 3, 4, and 5 provide the average tour length, gap to the optimal (provided by Concorde (Applegate et al., 2004)), and computational time, respectively. For space reasons, we only include results of eMAGIC(G) and eMAGIC(S) in Tables 1
and 2, and leave results of other versions in Table 4 in Appendix E.1. Moreover, we provide the result variances for our methods and evaluate them on TSP10,000 in resp. Appendices E.2 and E.3.
Tables 1 and 2 show that the computational times of exact solvers become prohibitive as the TSP size increases. Note that Gurobi is not able to solve TSP instances larger than 200 under a reasonable time budget. The classic heuristic methods are relatively fast, but their performances are not satisfactory. Among all learning methods, Fu et al. (2021) provide excellent results, but it is based on MonteCarlo tree search, which adapts to the instance to be solved. Without this adaptivity, our model with or without sampling provides competitive results. It can be better than Fu et al. (2021) up to TSP200. This somewhat suggests the limit of our approach, which trains on small instances and directly generalizes to large ones, without learning on the instance to be solved. We expect that our approach could be improved by training on slightly bigger instances (e.g., up to 100) or adding an adaptive component like in Fu et al. (2021).
Performance on Realistic TSP Table 3 compares the performances of a variety of learning-based TSP solvers on instances from TSPLIB. Each column of Table 3 represents the average gap to the optimal solution over the instances indicated by the corresponding rows. We treat the tested instances as small instances if their sizes are under 200 and large instances otherwise. And we extend the testing to larger instances with size up to 1002 and leave the performances of other models as empty for size ranging from 400 to 1002 since in their papers, the testings stop at instances with sizes around 400. We can observe from Table 3 that our model can perform much better than the other learning-based solvers not only for small problems but also large ones, which demonstrates the strong ability and the practical significance of eMAGIC to tackle realistic TSP problems. When the testing size increases, most models suffer from a relatively big increasing of the average gap while ours only increases less than 1% and remains in a good absolute value (3.40%) for even larger instances, which indicates a strong generalization ability of our model. Further details and experimental results are provided in Appendix I.
Ablation Study We demonstrate the strength of all the techniques we applied (including equivariance, stochastic CL, policy rollout baseline, combined local search, and RL) using an ablation study. We turn off each feature one at a time to see if the performance drops compared to the full version of eMAGIC(s1). For space reasons, details of the ablation study are deferred to Appendix G. Table 10 and 11 demonstrate that each technique plays a role in our model.

heading: 7 CONCLUSION
text: We presented a combination of novel techniques (notably, equivariance, combined local search, smoothed policy gradient) for designing an RL solver for TSP, which shows a good generalization capability. We demonstrated its effectiveness both on random and realistic instances, which shows that our model can reach state-of-the-art performance.
For future work, the approach can be further improved in various ways, e.g., extending it to the actor-critic scheme and exploiting invariances with the critic; or making it adaptive and learn on the instance to be solved. The approach could also be applied on other combinatorial optimization problems and other RL problems.

heading: A COMBINED LOCAL SEARCH
text: We first provide a brief introduction to TSP heuristics and then give the detailed presentation of each local search method we used, together with our whole combined local search algorithm.
Since TSP is an NP-hard problem (Papadimitriou, 1977), solving exactly large TSP instances is generally impractical. Hence, many heuristics have been proposed to solve TSP. Two important categories of heuristics are constructive heuristics and local search heuristics. On the one hand, constructive heuristics iteratively build a tour from scratch and therefore, do not rely on existing tours. As an example, insertion heuristics are constructive: they first choose a starting city randomly, then repeatedly insert an unvisited city into the partial tour that minimizes the increase of tour length until all the cities are included in the tour. The unvisited city can be selected in different ways leading to various versions of insertion heuristics: random insertion, nearest insertion, or farthest insertion (see Kool et al. (2019) for implementation details). Among them, farthest insertion usually yields the best results. On the other hand, local search heuristics try to improve a given complete tour by perturbing it. They are widely used as post optimization for TSP solvers, such as Deudon et al. (2018); Ma et al. (2019). One important class of local search heuristic is k-opt, which improves an existing complete tour σ by repeatedly performing the following operation: remove k edges of the current tour and reconnect the obtained subtours in order to decrease the tour length. For instance, one 2-opt operation would replace:
σ = ( σ(1), σ(2)..., σ(i), ..., σ(j), ..., σ(N) ) by
σ′ = ( σ(1), ..., σ(i), σ(j), σ(j − 1), ..., σ(i+ 1), σ(j + 1), ..., σ(N) ) where i < j < N if Lσ′(X) < Lσ(X). Based on k-opt, the LKH algorithm (Helsgaun, 2017) can often achieve nearly optimal solutions, but requires a long runtime.
Based on how many edges are removed and how the subtours are reconnected during one step, every local search heuristic can be seen as a special case of k-opt. Since the number of possible k-opt operations is O(Nk), 2-opt and 3-opt are usually preferred to efficiently search for fast improvements of existing tours, although they can get stuck in local optima. Hence, in our work, we use a combination of several local search methods, including local insertion heuristics, random 2-opt, search 2-opt, and search random 3-opt (see below for their definitions), to efficiently find potential 2-opt and 3-opt operations to improve the solution proposed by the RL solver. Instead of running each algorithm once for long enough, we run every local search shortly one by one and repeat this process for multiple times so that the combined local search are able to avoid more local optima. Intuitively, the rationale is that once a certain heuristic gets stuck in a local optimum, another can help get it out by trying a different operation. Next, we explain the above four heuristics, and then illustrate how we leverage them to form our combined local search.
A.1 RANDOM 2-OPT
For random 2-opt, we randomly try two edges for performing a 2-opt operation and repeat this process for α×Nβ times, where α, β > 0 are two hyperparameters controlling the strength of this heuristic as well as its runtime. Theoretically, random 2-opt can potentially cover all 2-opt operations and such a procedure makes random 2-opt much more flexible than trying all N(N − 1)/2 possible pairs of edges. In this paper, we set α = 0.5 and β = 1.5 for all experiments.

heading: A.2 LOCAL INSERTION HEURISTIC
text: For local insertion heuristic, inspired by the insertion heuristic, we iterate through all cities and find the best positions to insert them from their original locations. Let σ be the current tour and σt,t′ be the tour where we exchange the positions of σ(t) with σ(t′). Namely,
σt,t′ = ( σ(1), ..., σ(t′), σ(t), σ(t′ + 1), ..., σ(t− 1), σ(t+ 1), ..., σ(N) ) . (11)
for t′ 6= t − 1 and σt,t−1 = σ. The local insertion heuristic (see Algorithm 1) iterates over every t ∈ [N ]. For each t, we need to find t∗ such that:
t∗ = argmin t′ Lσt,t′ (X), (12)
where the definition of L is given in Equation (1). Then we replace σ by σt,t∗ . Theoretically, the local insertion heuristic is a special case of 3-opt where two of the removed edges cover a same node.
Algorithm 1 Local Insertion Heuristic 1: Input: A matrix of city coordinatesX = (xi)i∈[N ], current tour σ 2: Output: An improved tour σ 3: for t = 1 to N do 4: t∗ = argmint′ Lσt,t′ (X) 5: σ ← σt,t∗ 6: end for 7: return σ
A.3 SEARCH 2-OPT
For search 2-opt, we first iterate through all edges and for each first edge, we search for the best second edge to do the 2-opt. Let σ be the current tour and σ(t,t′) be the tour where we reverse the cities between σ(t) and σ(t′). Namely,
σ(t,t′) = ( σ(1), ..., σ(t− 1), σ(t′), σ(t′ − 1), ..., σ(t+ 1), σ(t), σ(t′ + 1), ..., σ(N) ) . (13)
where t < t′, and σ(t,t) = σ. Search 2-opt (see Algorithm 2) iterates over every t ∈ [N ]. For each t, we find t∗ such that:
t∗ = argmin t′≥t Lσ(t,t′)(X) (14)
Then we replace σ by σ(t,t∗). Search 2-opt potentially covers all possible 2-opt operations.
Algorithm 2 Search 2-opt 1: Input: A matrix of city coordinatesX = (xi)i∈[N ], current tour σ 2: Output: An improved tour σ 3: for t = 1 to N do 4: t∗ = argmint′≥t Lσ(t,t′)(X) 5: σ ← σ(t,t∗) 6: end for 7: return σ
A.4 SEARCH RANDOM 3-OPT
For search random 3-opt, the algorithm first randomly picks two edges. Then for these two randomly-picked edges, similar to search 2-opt, we search to find the best third edge to apply 3- opt. Let σ be the current tour, edges from σ(t1) to σ(t1 + 1) and σ(t2) to σ(t2 + 1) be the two randomly picked edges and σ∗(t1,t2,t3) be the optimal tour that differs from σ only in edges from σ(t1) to σ(t1+1), σ(t2) to σ(t2+1) and σ(t3) to σ(t3+1). For randomly picked t1 and t2, search random 3-opt (see Algorithm 3) iterates over every t3 ∈ [N ] and find t∗3 such that:
t∗3 = argmin t3 Lσ∗ (t1,t2,t3) (X) (15)
We repeat this process for α × Nβ times, where α, β > 0 are the same two hyperparameters as in random 2-opt. We set α = 0.5, β = 1.5 for all experiments. Search random 3-opt potentially covers all possible 3-opt operation.

heading: A.5 COMBINED LOCAL SEARCH
text: With these 4 different local search heuristics, our combined local search applies them one by one sequentially and repeat this for I times (see Algorithm 4), where I is a hyperparameter, which we set to I = 10 for all experiments in this paper.
Algorithm 3 Search random 3-opt 1: Input: A matrix of city coordinatesX = (xi)i∈[N ], current tour σ, hyperparameters α and β 2: Output: An improved tour σ 3: for iter= 1 to αNβ do 4: Randomly pick t1, t2 ∈ [N ] such that t1 6= t2 5: t∗3 = argmint3 Lσ∗(t1,t2,t3)(X) 6: σ ← σ∗(t1,t2,t∗3) 7: end for 8: return σ
Algorithm 4 Combined Local Search Algorithm 1: Input: A matrix of city coordinates X = (xi)i∈[N ], current tour σ, hyperparameters α, β and I for local search.
2: Output: An improved tour σ 3: for t = 1 to I do 4: σ ← apply local insertion heuristic(X , σ) 5: σ ← apply random 2-opt on σ for αNβ times 6: σ ← apply search 2-opt on σ 7: σ ← apply search random 3-opt on σ for αNβ times 8: end for 9: return σ

heading: B POLICY GRADIENT OF EMAGIC
text: As promised in Section 5 - Smoothed Policy Gradient, we elaborate on the detailed mathematical derivation for Equation (8).
The objective function J(θ) in Equation (3) can be approximated with the empirical mean of the total rewards using B trajectories sampled with policy πθ:
J(θ) = −E[Lσ(X)] ≈ −ÊB [ N∑ t=1 r (b) t ] = − 1 |B| |B|∑ b=1 N∑ t=1 r (b) t , (16)
where ÊB represents the empirical mean operator, Lσ(X) is the tour length of σ output by the RL policy and r(b)t is the t-th reward of the b-th trajectory. The policy gradient to optimize J(θ) can be estimated by:
∇θJ(θ) = −Eτ [( N∑ t=1 ∇θ log πθ(at|st) )( N∑ t=1 rt )]
≈ −ÊB [( N∑ t=1 ∇θ log πθ(a(b)t |s (b) t ) )( N∑ t=1 r (b) t )] (17)
However, Recall J(θ) is the standard objective used in most deep RL methods applied to TSP. Instead, we optimize J+(θ) = −E[Lσ+(X)] where Lσ+(X) is the tour length of σ after applying local search. This helps integrate better RL and local search by smoothing the value landscape and training an RL agent to output a tour that can be improved by local search. This new objective function can be rewritten:
J+(θ) = −Eσ∼πθ,σ+∼ρ(σ)[Lσ+(X)] = −Eσ∼πθ [Eσ+∼ρ(σ)[Lσ+(X)|σ]]
(18)
where ρ(σ) denotes the distribution over tours induced by the application of the stochastic local search on σ. Taking the gradient of this new objective:
∇θJ+(θ) = −Eτ [( N∑ t=1 ∇θ log πθ(at|st) ) Eσ+∼ρ(σ)[Lσ+(X)|σ] ]
≈ −ÊB [( N∑ t=1 ∇θ log πθ(a(b)t |s (b) t ) )( L σ (b) + (X(b)) )]
≈ − 1 |B| |B|∑ b=1 ( N∑ t=2 ∇θ log πθ(a(b)t |s (b) t ) )( L σ (b) + (X(b))− l(b) ) ,
(19)
where τ = (s1, a1, . . .) and σ is its associated tour. We simply approximate the conditional expectation over ρ(σ) by a sample. Therefore, our gradient estimate is an unbiased estimator of the gradient of our new objective J+(θ).
Using our policy rollout baseline introduces some bias to the estimation of the smoothed policy gradient. However, the variance reduction helps with achieving greater performance and smaller variance, as we observed in our experiments.

heading: C OVERALL ALGORITHM
text: We provide the pseudo-code of our overall training algorithm in Algorithm 5.
Algorithm 5 REINFORCE exploiting stochastic CL, equivariance, and smoothed policy gradient 1: Input: Total number of epochs E, training steps per epoch T , batch size B, hyperparameters α, β, γ and I for local search
2: Initialize θ 3: for e = 1 to E do 4: N ← Sample from p(e) according to Stochastic CL 5: for t = 1 to T do 6: ∀b ∈ {1, ..., B}X(b) ← Random TSP instance with N cities 7: ∀b ∈ {1, ..., B}σ(b) ← Apply πθ onX(b) after all the equivariant preprocessing steps 8: ∀b ∈ {1, ..., B}σ(b)+ ← Apply the combined local search on σ(b)
9: Use σ(b) and σ(b)+ to calculate∇θJ+(θ) according to Equation (8) 10: θ ← Update in the direction of∇θJ+(θ) 11: end for 12: end for

heading: D SETTINGS AND HYPERPARAMETERS
text: All our experiments are run on a computer with an Intel(R) Xeon(R) E5-2678 v3 CPU and a NVIDIA 1080Ti GPU. In consistency with previous work, all our randomly generated TSP instances are sampled in [0, 1]2 uniformly. During training, stochastic CL chooses a TSP size inR = {10, 11, . . . , 50} for each epoch e according to Equations (9) and (10), with σN set to be 3. We trained for 200 epochs in each experiment, with 1000 batches of 128 random TSP instances in each epoch. We set the learning rate to be 10−3 and the learning rate decay to be 0.96 in each experiment. In each experiment, we set α = 0.5, β = 1.5, γ = 0.25 and I = 10 for the parameters of our combined local search; As for random TSP testing and the ablation study, we test on TSP instances with size 20, 50, 100, 200, 500 and 1000 to evaluate the generalization capability of our model; as for realistic TSP, we test on TSP instances with sizes up to 1002 from the TSPLIB library. With respect to our model architecture, our MLP encoder has an input layer with dimension 2, two hidden layers with dimension 128 and 256, respectively, and an output layers with dimension 128; for the GNN encoder, we set H = 128 and nGNN = 3.

heading: E MORE EXPERIMENTS ON EMAGIC

heading: E.1 MORE VERSIONS OF EMAGIC
text: As promised in Section 6, we illustrate all versions of eMAGIC in this section:
E.2 VARIANCE ANALYSIS OF EMAGIC
As promised in Section 6, we provide the variance analysis of eMAGIC in this section. In our experiments, we repeated all our experiments (training + testing) with 3 random seeds. The variances are shown below:
The variances are quite low, showing that our method gives relatively good and stable results. Note the variances increase with the number of sampling (s1, s10 and S) for larger TSP instances since there is more room for improvement.

heading: E.3 EXPERIMENTS ON EXTREMELY LARGE TSP INSTANCES
text: In Table 7, we show the performances of eMAGIC on extremely large TSP instances (i.e., TSP 10,000) and the comparisons with a few methods that are able to generalize to TSP 10,000. For the hyperparameters of the combined local search, we use I = 2 and β = 1.4. We modify the
hyperparameters in this way because we find that for large TSP instances, doing too many iterations of local search is not efficient. The other experimental settings in this test are the same as Section 6. Also, we test Fu et al. (2021)’s method with a limited the time budget, denoted by AGCRN+Mlim: since we only reproduce their method successfully with CPU, we decrease their hyperparameter T (the MCTS will end no longer than T seconds) from 0.04n to 0.04n/1.8 ∗ (28/60) = 0.010n (1.8h and 28m are respectively the runtimes of their model and eMAGIC(G)). By this, we expect that their model and eMAGIC(G) will have a similar running time.
We can observe that our models outperform AM (Kool et al., 2019) to a large extent. Comparing to AGCRN+M (Fu et al., 2021) and LKH3 (Helsgaun, 2017), our methods run much faster without a big gap of performance. Plus, if we give AGCRN+M (Fu et al., 2021) a time budget comparable to our method (i.e., run AGCRN+M and eMAGIC(G) for the same amount of time), we can see that our algorithm outperforms AGCRN+M. Note that MCTS in AGCRN+M is written in C++. We could reduce our runtime if we had also written our code in C++ instead of Python. Therefore, our method offers a better trade-off in terms of performance vs runtime: it can generate relatively good results in much less time.

heading: F EVALUATION OF EQUIVARIANT PREPROCESSING METHODS

heading: F.1 COMPARISONS BETWEEN DIFFERENT SYMMETRIES USED DURING PREPROCESSING
text: We present the comparison results of applying rotation, translation and reflection. The comparisons are done by testing on random TSP instances followed the same setting in the experiment section. As in table 8, the rotation has the best performance on small and large TSP instances. By this, we choose rotation for our final algorithm.

heading: F.2 COMPARISONS BETWEEN ONE PREPROCESSING APPLICATION AND ITERATION PREPROCESSING APPLICATIONS
text: We present the comparison results of one preprocessing application and iteration preprocessing applications. The comparisons are done by testing on random TSP instances followed the same setting in the experiment section. As in table 9, the iteration preprocessing applications has better performance on small and large TSP instances. By this, we choose iteration preprocessing applications for our final algorithm.

heading: G ABLATION STUDIES
text: This section illustrates the details of Section 6 - Ablation Study.

heading: G.1 ABLATION STUDY ON KEY TRAINING TECHNIQUES
text: In this section, as promised in Section 6 - Ablation Study, we will illustrate the details of the ablation study regarding some key training techniques of eMAGIC. For equivariance, we remove all the equivariant procedures during training/testing (e.g., deleting visited cities, preprocessing, and using relative positions). For the policy rollout baseline, we replace it with the self-critic baseline, which is a greedy baseline implemented in Ma et al. (2019). For combined local search, we only apply it during testing to check if it can help improve the training process. For RL, we directly apply combined local search without performing any learning. Table 10 demonstrates that each technique plays a role in our model.

heading: G.2 ABLATION STUDY ON OTHER TRAINING TECHNIQUES
text: In this section, we first perform an ablation study of stochastic CL, meaning we fixed our TSP size to be 50 during training. Moreover, we perform an ablation study on each component of our equivalent model, which includes deleting the visited cities during training, equivariant preprocessing
operation, and using relative positions during training, as mentioned in Section 6 - Ablation Study. For the ablation study of these three components, we remove them during training and testing. Table 11 illustrates the results of our ablation study (including the full version for comparison), which demonstrate that each component is effective in our model.

heading: G.3 COMPARING PURE RL ALGORITHM WITH FULL ALGORITHM
text: In Section 6, when doing the ablation study for the combined local search, we do not apply it during training, but we still apply it during testing to check if it can help improve the training process. Here, we compare the pure RL algorithm (full algorithm without local search both in training and testing) and the full algorithm to demonstrate the strength of Combined Local Search:

heading: H GNN ENCODER DETAILS
text: Figure 2 is a detailed illustration of the GNN architecture introduced in Section 4.2. Please notice that Figure 2 could be regarded as a zoom-in version of the GNN part in Figure 1. The aggregation function used in GNN is represented by a neutral network followed by a ReLU function on each entry of the output.

heading: I REALIST TSP INSTANCES IN TSPLIB
text: As promised in Section 6, Tables 14 to 16 list the performances of various learning-based TSP solvers and heuristic approaches upon realistic TSP instances in TSPLIB. The bold numbers show the best performance among all the approaches. We can observe that most bold numbers are provided by eMAGIC, meaning our approach provides excellent results for TSPLIB. In addition, we also provide our results for the larger instances (with size up to 1002) from TSPLIB in Table 13. Table 13 only contains the experiments of our model since no other model can generalize to large size TSP instances with adequate performance.

Review summaries:
This paper introduces a deep RL approach combined with an equivariant model (to handle Euclidean symmetry) and local search heuristics (to improve a tour) to solve traveling salesman problems (TSP), in particular focusing on the generalizability of large-scale instances. The model consists of a graph neural network (GNN), a multi-layer perceptron (MLP), and an attention mechanism. In the training part, this model involves smoothed policy gradient, and stochastic curriculum learning to speed up the training and make the policy more generalizable. The experiments results show that the proposed approach significantly outperforms most learning-based solvers on large-scale randomly generated TSP and realistic TSP.

The paper proposes a series of techniques and a new model to learn to solve the TSP using RL. First, a number of preprocessing steps are presented in order to transform a TSP instance into a standardized form. Then, a new encoder-decoder architecture is proposed, based on GNN, MLP and attention modules. The authors propose a modification of the policy gradient algorithms for training by replacing the value of the policy-generated solution by its value after applying a local search heuristic, and use the original value as a baseline. Finally they use a curriculum learning approach to increase the difficulty of the instances seen during training, the difficulty here being represented by increasing the size of the instances, within a given interval. Experimental results are presented on standard synthetic and realistic TSP instances of size up to 1000 nodes.  

This paper proposes a novel combination of policy gradient and local search algorithms for the traveling salesman problems (TSP). The main idea is to apply local search to tours generated by policy rollout and compute policy gradients using the (potentially) improved tour from local search. In addition to the algorithmic contribution, this paper also presents a collection of preprocessing steps to ensure problem instance features are equivalent. Empirical studies are provided on random TSP instances as well as those in TSPLIB. Finally, ablation studies show the importance for each component of the proposed algorithm.

The paper employs equivariance properties and local search heuristics to derive a deep learning model capable of generalising to TSP instances of different sizes. The model constructs a tour, one city at a time, by learning a distribution (policy) over unvisited cities which is optimised via REINFORCE. The solution is further refined with a combination of different local search heuristics, and policy gradients are computed with respect to the final improved solution to smooth the loss landscape. The equivariance in TSP is mostly exploited during preprocessing steps—where rotation, translation and reflection operations are applied to the coordinates of each city—but is preserved by the model architecture which consists of a graph neural network and attention mechanisms. The model is trained with stochastic curriculum learning only on small TSP instances ($\leq 50$ cities) but is shown to have competitive performance on problems with up to 1000 cities.