Unsupervised Learning for Solving the Travelling Salesman Problem

Given a TSP instance, let $\mathbf{D}_{i,j}$ denote the Euclidean distance between city $i$ and city $j$. $\mathbf{D}\in\mathbb{R}^{n\times n}$ is the distance matrix. We first build adjacency matrix $\mathbf{W}\in\mathbb{R}^{n\times n}$ with $\mathbf{W}_{i,j}=e^{-\mathbf{D}_{i,j}/\tau}$ and node feature $\mathbf{F}\in\mathbb{R}^{n\times 2}$ based on the input coordinates, where $\mathbf{F}_{i}=(x_{i},y_{i})$ and $\tau$ is the temperature. The node feature matrix $\mathbf{F}$ and the weight matrix $\mathbf{W}$ are then fed into a GNN to generate a soft indicator matrix $\mathbb{T}\in\mathbb{R}^{n\times n}$.

In our model, we use Scattering Attention GNN (SAG), SAG has both low-pass and band-pass filters and can build adaptive representations by learning node-wise weights for combining multiple different channels in the network using attention-based architecture. Recent studies show that SAG can output expressive representations for graph combinatorial problems such as maximum clique and remain lightweight Min et al. (2022).

Let $\mathcal{S}\in\mathbb{R}^{n\times n}$ denote the output of SAG, we first apply a column-wise Softmax activation to the GNN’s output and we can summarize this operation in matrix notation as $\mathbb{T}_{i,j}={e^{\mathcal{S}_{i,j}}}/{\sum_{k=1}^{n}e^{\mathcal{S}_{k,j}}}$. This ensures that each element in $\mathbb{T}$ is greater than zero and the summation of each column is 1. We then use $\mathbb{T}$ to build a heat map $\mathcal{H}$, where $\mathcal{H}\in\mathbb{R}^{n\times n}$.

In our model, we use $\mathcal{H}$ to estimate the probability of each edge belonging to the optimal solution and use $\mathbb{T}$ to build a surrogate loss of the Hamiltonian Cycle constraint. This will allow us to build a non-smooth heat map $\mathcal{H}$ and improve the performance of the local search.

Before building the unsupervised loss, let’s recall the definition of TSP. The objective of TSP is to identify the shortest Hamiltonian Cycle of a graph. Therefore, the unsupervised surrogate loss should act as a proxy for two requirements: the Hamiltonian Cycle constraint and the shortest path constraint. However, designing a surrogate loss for the Hamiltonian Cycle constraint can be challenging, particularly when working with a heat map $\mathcal{H}$. To address this, we introduce the $\mathbb{T}\rightarrow\mathcal{H}$ transformation, which enables the model to implicitly encode the Hamiltonian Cycle constraint.

To better understand the $\mathbb{T}\rightarrow\mathcal{H}$ transformation, we show a binary instance in Figure 1. Figure 1 illustrates a soft indicator matrix $\mathbb{T}$, its heat map $\mathcal{H}$ following the transformation $\mathbb{T}\rightarrow\mathcal{H}$, and their corresponding routes. When we directly use the soft indicator matrix $\mathbb{T}$ as the heat map. It can result in loops (parallel edges) between cities, such as (2,3) and (4,5) in Figure 1 (middle). After we apply the $\mathbb{T}\rightarrow\mathcal{H}$ transformation, the corresponding heat $\mathcal{H}$ is a Hamiltonian Cycle, as shown in the right part in Figure 1.

We build the heat map $\mathcal{H}$ based on $\mathbb{T}$. As mentioned, $\mathcal{H}_{i,j}$ is the probability for edge ($i$,$j$) to belong to the optimal TSP solution. We define $\mathcal{H}$ as:

where $p_{t}\in\mathbb{R}^{n\times 1}$ is the $t_{th}$ column of $\mathbb{T}$, $\mathbb{T}=[p_{1}|p_{2}|...|p_{n}]$. As shown in Figure 1, the first row in $\mathcal{H}$ is the probability distribution of directed edges start from city $1$, and since the third element is the only non-zero one in the first row, we then add directed edge $1\rightarrow 3$ to our TSP solution. Similarly, the first column in $\mathcal{H}$ can be regarded as the probability distribution of directed edges which end in city $1$. Ideally, given a graph $\mathcal{G}$ with $n$ nodes, we want to build a soft indicator matrix where each row and column are assigned with one value 1 (True) and $n-1$ values 0 (False), so that the heat map will only contain one valid solution. In practice, we will build a soft indicator matrix $\mathbb{T}$ whose heat map $\mathcal{H}$ assigns large probabilities to the edges in the TSP solution and small probabilities to the other edges.

Overall, the $\mathbb{T}\rightarrow\mathcal{H}$ transformation in Equation 1 enables us to build a proxy for the Hamiltonian Cycle constraint. We further prove that $\mathcal{H}$ represents one Hamiltonian Cycle when each row and column in $\mathbb{T}$ have one value 1 (True) and $n-1$ value 0 (False). We refer to the proof in appendix D.

In order to generate such an expressive soft indicator matrix $\mathbb{T}$, we minimize the following objective function:

The first term in $\mathcal{L}$ encourages the summation of each row in $\mathbb{T}$ to be close to 1. As mentioned, we normalize each column of $\mathbb{T}$ using Softmax activation. So when the first term is minimized to zero, each row and column in $\mathbb{T}$ are normalized (doubly stochastic). The second term penalizes the weight on the main diagonal of $\mathcal{H}$, this discourages self-loops in TSP solutions. The third term can be regarded as the expected TSP length of the heat map $\mathcal{H}$, where $\mathbf{D}_{i,j}$ is the distance between city $i$ and $j$. As mentioned, since $\mathcal{H}$ corresponds to a Hamiltonian Cycle given an ideal soft indicator matrix with one value $1$ (True) and $n-1$ value $0$ (False) in each row and column. Then the minimum value of $\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{D}_{i,j}\mathcal{H}_{i,j}$ is the shortest Hamiltonian Cycle on the graph, which is the optimal solution of TSP. To summarize, we build a loss function which contains both the shortest and the Hamiltonian Cycle constraints. The shortest constraint is realized by minimizing $\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{D}_{i,j}\mathcal{H}_{i,j}.$ For the Hamiltonian Cycle constraint, instead of writing it in a Lagrangian relaxation style penalty, we use a GNN which encourages a non-smooth representation, along with the doubly stochastic penalty, and the $\mathbb{T}\rightarrow\mathcal{H}$ transformation.

Given a heat map $\mathcal{H}$, we consider $M$ largest elements in each row (without diagonal elements) and set other $n-M$ elements as 0. Let $\tilde{H}$ denote the new heat map, we then symmetrize the new heat map by $\mathcal{H}^{\prime}=\tilde{H}+\tilde{H}^{T}$. Let $\mathbf{E}_{ij}\in\{0,1\}$ denote whether an undirected edge ($i,j$) is in our prediction or not. Without loss of generality, we can assume $0<i<j\leq n$ and define $\mathbf{E}_{ij}$ as :

Let $\Pi$ denote the set of undirected edges $(i,j)$ with $\mathbf{E}_{ij}=1$. Ideally, we would build a prediction edge set $\Pi$ with a small $M$ value, and $\Pi$ can cover all the ground truth edges so that we can reduce search space size from $n(n-1)/2$ to $|\Pi|$. In practice, we aim to let $\Pi$ cover as many ground truth edges as possible and use $\mathcal{H}^{\prime}$ to guide the local search process.

We employ the best-first local search guided by the heat map to generate the final solution. Best-first search is a heuristic search that explores the search space by expanding the most promising node selected w.r.t. an evaluation function $f(node)$. In our framework, each node of the search tree is a complete TSP solution. For the initialization of one search tree, we randomly generate a valid TSP solution and improve it using the 2-opt heuristic until no better solution is found. The expand action of the search node refers to Fu et al. (2021) and is based on the widely used k-opt heuristic Croes (1958), where it replaces $k$ old edges (in the current solution) with $k$ new edges, i.e., transforms the old solution to a new solution. More formally, we use a series of cities, $u_{1},v_{1},u_{2},\ldots,u_{k},v_{k+1}$, to represent an action, where $v_{k+1}=u_{1}$ to ensure it is a valid solution. All the edges $(u_{i},v_{i})\ (1\leq i\leq k)$ are removed from the tour and $(v_{i},u_{i+1})(1\leq i\leq k)$ are added to the tour. Note that once we know $u_{i}$, $v_{i}$ is deterministically decided. $u_{1}$ is randomly selected for each expansion, and $v_{1}$ is decided subsequently. Then we select $u_{i+1}(i\geq 1)$ as follows: (1) if $u_{i+1}=u_{1}$, i.e., forms a new TSP tour, leads to an improved solution then we set $u_{i+1}=u_{1}$ and have a candidate solution. (2) if $i\geq K$, then we will discard this action and start a new expand action, where $K$ is a hyper-parameter which controls the maximal edges we can remove in one action. (3) otherwise, we select $u_{i+1}$ based on the heat map stochastically. We use $N_{u,v}$ to denote the times the edge $(u,v)$ is selected during the entire search procedure. The likelihood of selecting the edge $(u,v)$ is denoted by $L_{u,v}\ =\ \mathcal{H}^{\prime}_{u,v}+\alpha\sqrt{\frac{\log(S+1)}{N_{u,v}+1}}$, where $\alpha$ is a hyper-parameter and $S$ is the local search’s total number of expand actions. The first term encourages the algorithm to select the edge with a high heat map value, while the second term diversifies the selected edges. Moreover, when selecting the city $v$ given $u$, we only consider the cities from the candidate set of $v$. This candidate set consists of cities with the top $M$ heat map value or the nearest $M$ cities

Among all the possible new solutions, we use the tour’s length as the evaluation function $f$, i.e., we select the solution of the shortest tour length as the next search node. For each search node, we try at most $T$ expand actions. From these $T$ expand actions, if no improved solution is found, we randomly generate a new initial solution and start another round of best-first local search.

We borrow the idea of the backpropagation used in Monte Carlo Tree Search (MCTS). We use $s$ to denote the current search node, $s^{\prime}$ to denote the next search node ($s^{\prime}$ has to improve $s$), and $L(s)$ to represent the tour length of node $s$. The heat map $\mathcal{H}^{\prime}$ is updated as:

where $\beta$ is a search parameters and $(v_{i},u_{i+1})(1\leq i\leq k)$ is the actions used to transform $s$ to $s^{\prime}$. We raise the importance of the edges that lead to a better solution. If we cannot find an improved solution for the current node, then no update is executed for the heat map.

Randomness is shown to be very powerful in the local search community Gomes et al. (1998); Bresina (1996); Gomes et al. (2000). Our local search procedure also employs it to improve performance. When we stop the current round of local search by not finding an improved solution within $T$ expand actions and switch to a new best-first local search with a new initial solution, we randomly modify the parameter $K$. A larger $K$ value results in more time searching from one initial solution. The intuition is that sometimes we want more initial solutions while sometimes we want to search deeper (replace more edges in k-opt) for a specific solution. Besides that, we also randomly decide how we construct the candidate set (based on the heat map or pairwise distance) for each city before the new round of local search.

Our dataset contains 2,000 samples for training and 1,000 samples for validation. We use the same test dataset in Fu et al. (2021). The test dataset contains $10,000$ 2D-Euclidean TSP instances for $n=20,50,100$, and 128 instances for $n=200,500,1,000$. We train our models on TSP instances with 20, 50, 100, 200, 500, and 1,000 vertices. We then build the corresponding heat maps based on these trained models.

Table 1, Table 2 and Table 3 present model’s performance on TSP 20, 50, 100, 200, 500 and 1,000. The first three lines in the tables summarize the performance of two exact solvers (Concorde and Gurobi) and LKH3 heuristic Helsgaun (2017). The learning-based methods can be divided into RL sub-category and SL sub-category. Greedy decoding (G), Sampling (S), Beam Search (BS), and Monte Carlo Tree Search are the decoding schemes used in RL/SL. The 2-OPT is a greedy local search heuristic.

We compare our model with existing solvers as well as different learning-based algorithms. The performance of our method is averaged of four runs with different random seeds. The running time for our method is divided into two parts: the inference time (building the heat map $\mathcal{H}$) and the search time (running search algorithm).

On small instances, our results match the ground-truth solutions and generate average gaps of -0.00009%, -0.002% and -0.00011% respectively on instances with $n=20,50,100$, where the negative values are the results of the rounding problem. The total runtime of our method remains competitive w.r.t. all other learning baselines. On larger instances with $n$ = $200,500$ and $1,000$, we notice that traditional solvers (Concorde, Gurobi) fail to generate the optimal solutions within reasonable time when the size of problems grows. For RL/SL baselines, they generate results far away from ideal solutions, particularly for cases with $n=1,000$. Our UTSP method is able to obtain 0.0918%, 0.8394% and 1.1770% on TSP $200,500$ and $1,000$, respectively. We remark that UTSP outperforms the existing learning baselines on larger instances (TSP 200, 500, 1000) ¹¹1On TSP 500, when we increase the time budget for Search, we achieve 0.42% in 1.37m + 8.32m.. More discussion between Fu et al. (2021) and UTSP can be found in Appendix C.

Our model takes less training time than RL/SL method because we require very few training instances. Taking TSP 100 as an example, RL/SL needs 1 million training instances, and the total training time can take one day using a NVIDIA V100 GPU, while our method only takes about 30 minutes with 2,000 training instances. The training data size does not increase w.r.t. TSP size. Our training data consists of 2,000 instances for TSP 200, 500 and 1,000. At the same time, the UTSP model also remains very lightweight. On TSP 100, we use a 2-layer SAG with 64 hidden units and the model consists of 44,392 trainable parameters. In contrast, RL method in Kool et al. (2019) takes approximately 700,000 parameters and the SL method in Joshi et al. (2022) takes approximately 350,000 parameters.

Our UL method generalizes well to unseen examples without requiring a large number of training samples. This is because the loss function in Equation 2 is fully differentiable w.r.t. the parameters in SAG and we are able to train the model in an end-to-end fashion.

In other words, given a heat map $\mathcal{H}$, the model learns to assign large weights to more promising edges and small weights to less promising ones through backpropagation without any prior knowledge of the ground truth or any exploration step. However, when using SL, the model learns from the TSP solutions, which fails when multiple solutions exist or the solutions are not optimal Li et al. (2018). While for RL, the model often encounters an exploration dilemma and is not guaranteed to converge Bengio et al. (2021)Joshi et al. (2019b). Overall, UTSP requires fewer training samples and has better generalization compared to SL/RL models.

We aim to generate a non-smooth soft indicator matrix $\mathbb{T}$ and build an expressive heat map $\mathcal{H}$ to guide the search algorithm.

However, most GNNs aggregate information from adjacent nodes and these aggregation steps usually consist of local averaging operations, which can be interpreted as a low-pass filter and causes the oversmoothing problem Wenkel et al. (2022). The low-pass model generates a smooth soft indicator matrix $\mathbb{T}$, which finally makes the elements in $\mathcal{H}$ become indistinguishable. So it becomes difficult to discriminate whether the edges belong to the optimal solution or not. In our model, we assume all nodes in the graph are connected, so every node has $n-1$ connected to neighbouring nodes. This means every node receives messages from all other nodes and we have a global averaging operation over the graph, this can lead to severe oversmoothing issue.

To avoid oversmoothing, one solution is to use shallow GNNs. However, this would result in narrow receptive fields and create the problem of underreaching Barceló et al. (2020). Our model uses SAG because this scattering-based method helps overcome the oversmoothing problem by combining band-pass wavelet filters with GCN-type filters Min et al. (2022). Figure 2 illustrates the training loss on TSP $100$ and the differences between our SAG model and the graph convolutional network (GCN) Kipf and Welling (2016), where GCN only performs low-pass filtering on graph signals Nt and Maehara (2019). When using GCN, the training loss decreases slowly, and the validation loss reaches a plateau after we train the model for 20 epochs. This is because the low-pass model generates a smooth $\mathbb{T}$. Such a smooth $\mathbb{T}$ results in an indistinguishable $\mathcal{H}$, detrimentally impacting the training process. Instead, we observe lower training and validation loss when using SAG; this suggests that SAG generates a more expressive representation which facilitates the training process.

Figure 3 illustrates the generated heat maps using GCN and SAG on a TSP 100 instance, we choose this instance from the validation set randomly. When using the GCN, due to the oversmoothing problem, the model generates a smooth representation and $\mathcal{H}$ becomes indistinguishable. The elements in $\mathcal{H}$ have a small variance and most of them are $\sim$ 0.01. Instead, the SAG generates a discriminative representation and the elements in the heat map have a larger variance.

Here, we train both GCN and SAG with the same loss function. So the differences illustrated in Figure 3 are the direct result of overcoming the oversmoothing problem.