This article demonstrates the effectiveness of employing a deep learning model in an optimization pipeline. Specifically, in a generic exact algorithm for a NP problem, multiple heuristic criteria are usually used to guide the search of the optimum within the set of all feasible solutions. In this context, neural networks can be leveraged to rapidly acquire valuable information, enabling the identification of a more expedient path in this vast space. So, after the explanation of the tackled traveling salesman problem, the implemented branch and bound for its classical resolution is described. This algorithm is then compared with its hybrid version termed "graph convolutional branch and bound" that integrates the previous branch and bound with a graph convolutional neural network. The empirical results obtained highlight the efficacy of this approach, leading to conclusive findings and suggesting potential directions for future research.

§ INTRODUCTION

The well-known Traveling Salesman Problem (TSP) calls for finding a minimum cost tour
among $n$ cities. In the most common formalization, an undirected
graph $G(V,E)$ is given, with each vertex in $V=\{1,2,\dots,n\}$ representing a city
and each edge $e=\{i,j\}\in E$ (carrying a weight $w_e$) representing the cost of moving from city $i$ to $j$.
The problem requires to determine a Hamiltonian cycle $H^\star$ passing through
all vertices of $G$ exactly once so that the total cost $w(H^\star)=\sum_{e\in H^\star} w_e$ is
as small as possible.

The problem has become a classical one in Combinatorial Optimization (CO), and an enormous amount of research has been devoted to its solution. From a theoretical point of view, the problem is NP-complete. However, in practice it is among the best-solved ones, with exact algorithms that

are able to deliver optimal solutions on instances with thousands of cities; see <cit.> for a comprehensive study.
Excellent heuristic algorithms, which deliver high-quality solutions in short computation times are also available, for instance the Kerinighan-Lin local search <cit.>.

Despite this, the amount of active research on the TSP is still growing, on one hand because of its wide applicability (from network design to logistics, with more and more demanding large-scale applications).
On the other hand, precisely because of the availability of such excellent algorithms, it has become a favourite test-subject for almost every new optimization technique.

In the so-called metric variant (metric TSP, m-TSP) of the problem the graph is complete, with edge weights obeying the triangular
inequality $w_{\{i,j\}}+w_{\{j,k\}}\geq w_{\{i,k\}}$ for each triple $i,j,k\in V$. The m-TSP is still NP-complete
but approximable within a constant relative error <cit.>.
We consider the geometric TSP (g-TSP) where each city/graph vertex is associated with a point in the bi-dimensional euclidean plane.

In this paper we investigate the possibility of enhancing exact and heuristic procedures for the g-TSP via the introduction of Machine-Learning (ML) techniques.
Machine Learning deals with the problem of building systems that are able to improve their performances on a particular task given some previous experience <cit.>.

In the last decade, ML has found application across a wide range of problems, gaining global prominence, particularly propelled by the remarkable successes attained in Deep Learning (DL) <cit.>.

Recently, there have been efforts to integrate traditional optimization algorithms with ML models. For example, one of the first attempts that tries to integrate machine learning in an optimization algorithm is <cit.> that learns a ranking function to approximate the decisions made by the strong branching in a branch and bound solver. A comprehensive survey on the application of machine learning to combinatorial optimization problems is presented by <cit.>. While DL is not suitable for exact solutions to CO problems, it can enhance the performance of classical exact or heuristic algorithms.

In this work, we exploit a deep learning system to predict the probability that edges of a graph belong to a g-TSP solution.

Following this direction, in this paper we introduce a novel hybrid exact algorithm, called Graph Convolutional Branch and Bound (GCBB), for the g-TSP. GCBB integrates a Graph Convolutional neural Network (GCN) with a traditional Branch and Bound (BB) algorithm.
The computational overhead of using a deep neural network is justified by the effective utilization of extracted information, resulting in faster computations and in the exploration of fewer branch and bound nodes than the traditional approach. The operational workflows of the conventional branch and bound and its enhanced hybrid neural counterpart are illustrated in Figure <ref>.

Hence in the authors' view, the contribution of this paper is twofold.

* We show that there is room for incorporating information of the DL model
into an exact algorithm, enhancing its overall performances. We incorporate
knowledge from a trained DL model into a simple branch and bound for the TSP,
allowing to add some heuristic guidance to the exploration of the search tree,
and obtaining a significant improvement in performances.

* Training a DL model is computationally expensive, and often the resulting
model is specifically tailored to a given size of the instances.
We show that, with some simple manipulation, the heuristic information from the DL
model can still be profitably used without resorting to training a new model.

The first attempt to use neural networks for the TSP dates back to a Hopfield network <cit.>, while the pioneer paper that tried to use modern DL to tackle TSP and other combinatorial optimization problems is <cit.>. For problems like the TSP that concern graphs the use of a graph neural network capable of learning topological information of the involved instance has shown to be promising. An example of these are the Graph Attention Network  <cit.> and the Graph Convolutional Network (GCN) <cit.>. Most of these cited works focus on entirely replacing an algorithm with a DL model, in contrast with the hybrid algorithm developed in this work.

The remainder of this paper is structured as follows. A classical branch and bound algorithm for the TSP is described in <ref>.
Then, its improved version which relies on a trained GCN is described in <ref> and <ref>. [Code, data and results are available at the following GitHub repository: <https://github.com/LorenzoSciandra/GraphConvolutionalBranchandBound>.]
[sec:Experiments]Section <ref>sec:Experiments is designated for detailing all the conducted experiments, while the analysis of computational results, comparing the two algorithms, is provided in <ref>.
In <ref> we offer conclusions with possible directions for future research.

§ BRANCH AND BOUND FOR THE GEOMETRIC TSP

The starting point for solving the g-TSP, which is enhanced in <ref>, is the classical branch and bound scheme proposed by Held and Karp <cit.>, with some improvements from more recent works <cit.>.

The g-TSP is formalized by means of a complete graph $G=(V,E)$ where $V$ is the set of $n$ vertices and $E$ is the set of the $m = \frac{n\cdot(n-1)}{2}$ edges. We denote with $\delta(i)=\{\{i,j\} \in E \mid j\in V\}$ the set of edges incident on vertex $i$, and with $E(S)=\{\{i,j\}\in E\mid i,j\in S \subset V\}$  the set of edges connecting pairs of vertices in a subset $S$ of $V$. The TSP model is the following:

\begin{align}
\min \:f(\boldsymbol x) = \min & \sum_{e\in E} w_ex_e\label{tsp:obj} \\
\text{subject to}\ & \sum_{e\in\delta(i)} x_e=2 && i\in V\setminus\{n\}\label{tsp:degree}\\
& \sum_{e\in E}x_e=n \label{tsp:num-edges}\\
& \sum_{e\in E(S)} x_e \leq |S| - 1  && \forall\: S \subset V \setminus\{1\}\label{tsp:subtour} \\
&x_e \in \{0,1\} &&\forall e\in E\label{tsp:vars}
\end{align}

Binary decision variables $\boldsymbol x=(x_e\colon e\in E)$ are defined so that $x_e=1$ iff edge $e$ is selected to be in the tour; constraints (<ref>)–(<ref>)
force every vertex to have exactly two selected incident edges.

Constraints (<ref>) are the so-called sub-tour elimination
constraints forbidding non-connected solutions.

Following <cit.>, we use a lower bound based on a Lagrangian relaxation combined with an iterative sub-gradient improvement, that achieves the optimum of the original g-TSP problem through branching iterations. Upon choosing a starting node, e.g. $1$, for computing the Lagrangian relaxation, degree constraints (<ref>) for $i\neq 1$ are
relaxed into the objective function via multipliers
$\boldsymbol{\pi}=(\pi_i:i\in V\setminus\{1\})$. The relaxed
problem is then reduced to the following:

\begin{align}
L(\boldsymbol\pi) = \min_{\boldsymbol x} f(\boldsymbol x;\boldsymbol\pi)= \min_x &
\sum_{e\in E} w_ex_e+\sum_{i\in V\setminus\{1\}}\pi_i\left(\sum_{i\in\delta(i)}x_e-2\right)
\tag{\ref{tsp:obj}'}\label{1-tree:obj} \\
\text{subject to}\ & \sum_{e\in\delta(1)} x_e=2\tag{\ref{tsp:degree}'}\label{1-tree:degree1}\\
& \sum_{e\in E}x_e=n \tag{\ref{tsp:num-edges}'}\label{1-tree:num_edges}\\
& \sum_{e\in E(S)} x_e \leq |S| - 1  && \mathllap{\forall\: S \subset V\setminus\{1\}}
\tag{\ref{tsp:subtour}'}\label{1-tree:subtour} \\
&x_e \in \{0,1\}&&\forall e\in E\tag{\ref{tsp:vars}'}\label{1-tree:vars}
\end{align}

The integer program (<ref>)–(<ref>), is solved
by computing a so-called minimum-cost 1-tree <cit.> on graph $G$ with reduced edge costs $w_e'=w_e+\pi_i+\pi_j$ (with $e=\{i,j\}$). This is done quickly with well-known polynomial

In order to determine the best possible multipliers that maximize the lower bound $L(\boldsymbol\pi),$
the vector $\boldsymbol\pi$ is iteratively adjusted by sub-gradient optimization accordingly with the procedure proposed in  <cit.> and  <cit.>.

All the previous procedures, with others that we now explain, are collected in the branch and bound algorithm presented in <ref> that we propose for the resolution of the g-TSP. It is worth noting that the algorithm takes as input a symmetric matrix $\mathbf{P}$ of dimension $|V| \times |V|$. Its presence is instrumental to the extension we introduce in <ref>. For the time being we assume that $\textbf{P}=0$, with the understanding that in so doing the presented algorithm reverts to a standard resolution of the TSP via a branch and bound approach.

Generalized Branch and Bound – if $\mathbf{P} = \mathbf{0}$, the algorithm reverts to the classical branch and bound algorithm. Otherwise the algorithm is guided by the heuristic information in $\mathbf{P}$.
Instance $G=(V,E)$, matrix of edges probabilities $\mathbf{P}$, and the maximum time $\texttt{time\_limit}$.
The optimal tour $T^\star$.
$T = \emptyset$*1-tree solution
$c = \texttt{OPEN}$*subroblem class
$\mathbf{v} = \texttt{starting\_vertex}(G,\mathbf{P})$
$S = (G, T, c, \mathbf{v})$ *the original problem
$S^\star = \texttt{first\_solution}(G,\mathbf{P})$  *the best solution found
$L=\{S\}$*ordered list of all subproblems
$L \neq \emptyset \: \land \: \texttt{curr\_time} < \texttt{time\_limit} $
$S = \texttt{pop(L)}$
$S^\star = \texttt{bound}(S, S^\star, \mathbf{P})$ *possibly close the node
$S.c == \texttt{OPEN}$
$\texttt{variable\_fixing}(S)$*possibly close the node
$\texttt{infer\_constraints}(S)$*possibly close the node
$S.c == \texttt{OPEN}$
$L_S = \texttt{branch}(S, \mathbf{P})$ *the new subproblems
$L = \texttt{add}(L, L_S, \mathbf{P})$

After selecting the starting vertex at <ref> by maximizing the lower bound to minimize the initial gap with the optimum, a multi-start nearest-neighbor heuristic (see <cit.>) is applied on <ref> at the root node to generate an initial feasible solution.
The repeated resolution of the relaxed problem is performed at <ref>, and, if the branch and bound node is not closed, some edge-fixing techniques from <cit.> (<ref> and <ref>)  are applied in order to further restrict whenever possible the number of edges that are candidate for the branching operation.
The branching scheme adopted in <ref> is (assuming $\textbf{P} = 0$) the one proposed in <cit.>, where an edge-scoring mechanism is coupled with the sub-gradient optimization procedure. The edge $e\in E$ with the “best Sthluer score” is selected for a binary branching operation: on one branch the edge $e$ is forced to belong to the
tour (setting $x_e=1$), in the other branch $e$ is forbidden

§ GRAPH CONVOLUTIONAL BRANCH AND BOUND

Since the prohibitive complexity of the g-TSP stems from the exponential number of required 1-trees to be explored, a key challenge is to devise a criterion that curbs this combinatorial explosion. To address this challenge, in this section, we introduce a score for choosing among equally promising 1-trees, i.e. when the difference in value of their lower bower bounds on the optimum is less than a reasonable significance threshold. We evaluate such score using the output of the GCN, which enables us to speed up the classical branch and bound subprocedures defining the new Graph Convolutional Branch and Bound that we present in <ref>.

In our approach, a graph neural network is trained in a supervised manner by taking as input the coordinates of all the vertices of the graph to produce in output the symmetric probability matrix $\mathbf{P}$, whose $p_{ij}$ entries are the probabilities that the edges $e_{ij}$ belong to an optimal tour. For this purpose, the GCN uses node and edge embeddings to convey topological information about the instance graph. Notably, the edge embeddings act as a gating mechanism, modulating the information conveyed by neighboring node embeddings. This process resembles a dense attention map (<cit.>), with the modulated information being added to the embedding of the previous node layer. Concurrently, the edge embeddings are derived by summing the embeddings of the previous layer, with a non-linear transformation of the node embeddings that the edge connects.

In all our experiments we leveraged the GCN
models trained by <cit.>.[Pre-trained pytorch models can be downloaded from <https://github.com/chaitjo/graph-convnet-tsp>.]

§.§ 1-Tree's optimality score

In order to propose a criterion for choosing among 1-trees

of a given instance graph $G$ we define the probability space $(\Omega, \mathcal{F}, \mathbb{P})$. Here $\Omega$ is the set of edges $e_{ij}$ connecting every pair of vertices, $\mathcal{F}$ is the $\sigma$- algebra that contains all subsets of $\Omega$ and $\mathbb{P}$ is the corresponding probability measure.
On this probability space we introduce the subset $R$ that includes all the edges that appear at least once in an optimal tour and we call $\bar{P}$ the matrix with elements $\bar{p}_{ij}=\mathbb{P}(e_{ij}\in R)$, for all $i$,$j$.
Since we are not able to compute $\bar{P}$, we approximate it using $\mathbf{P}$, which is fully characterized by the output of the GCN. However, since it is derived from the neural network and is contingent upon its training, although all optimal tours are theoretically equivalent, the GCN may exhibit a bias towards specific tours, deviating from a uniform distribution over them.

Let us now denote with $I_{ij}$ the indicator random variable taking the value $1$ when edge $e_{ij}$ is present in an optimal tour according to the neural network. Given the discussion above, we assume that $I_{ij}$ is equal to $1$ with probability $p_{ij}$. Following this line of reasoning, since the 1-tree is a relaxation of the tour structure, and the set of all 1-trees is a super-set of all possible tours, the previously defined indicator variables can be utilized to evaluate the proximity of a 1-tree $T$ to an optimal tour. To do that we introduce a new random variable $O_T = \sum_{e_{ij} \in T} I_{ij}$, termed the optimality score of the 1-tree $T$. The resulting value of $O_T$ represents the number of edges in the 1-tree $T$ that belong to at least one optimal tour, and its expected value can be computed as

\mathbb{E}[O_T] = \mathbb{E}\left[\sum_{e_{ij} \in T} I_{ij}\right] = \sum_{e_{ij} \in T} \mathbb{E}[I_{ij}] = \sum_{e_{ij} \in T} p_{ij} \in [0, n]

The criterion that we propose for choosing among different 1-trees $T_i$ is the maximization of the expected number of optimal edges that the 1-trees contain

\max_{T_{i}} \:\: \mathbb{E}[O_{T_i}] = \max_{T_{i}} \sum_{e_{ij} \in T} p_{ij}

and we use it when other optimization rules are unable to distinguish among them.

In the following sections we discuss the effect of this criterion in <ref> when the matrix $\mathbf{P}\neq \mathbf{0}$ is derived from the output of the neural network.

§.§ Starting vertex selection

Upon assigning the optimality score to each 1-tree, the selection of the starting vertex in the GCBB, <ref> in <ref>, is determined by choosing the one that attains the maximum lower bound or the one with the highest expected value $\mathbb{E}[O_T]$ when a unique maximum cannot be found. So, unlike the classical branch and bound method, which solely focuses on maximizing the lower bound to partially close the gap with the upper bound, in its neural version the information derived from the GCN is employed to break ties.

§.§ Probabilistic Nearest Neighbor

Leveraging the probabilistic information of the edges now available, we advocate utilizing, at <ref> in <ref>, an approach we term the Probabilistic Nearest Neighbor (PNN) algorithm for determining the initial feasible solution. The PNN starts from a specific vertex, and, at each step, chooses the edge with the maximum probability that leads to an unvisited vertex. This process is iteratively applied with all vertices as starting points. The minimum tour among all iterations is then compared with the one obtained through the multi-start Nearest Neighbor approach described in <ref>. Ultimately, only the best tour reached is regarded as the initial feasible solution.

§.§ Neural Branch

Among edges that attain the same best value of the Shutler's score, the algorithm selects the branching edge (<ref> in <ref>) with the highest probability $p_{ij}$.
The rationale behind this is that if the GCN assigns a high probability to an edge because it deems it likely to be part of an optimal tour, then, among peers with the same Shutler's score, it is sensible to prioritize exploration of this edge first.

§.§ Node Selection

In a branch and bound, when a subproblem with a value lower than the current upper bound generates its children through the branching step, these children are added into a list of open decision nodes (<ref> in <ref>). In the classic version of the algorithm. i.e. when $\mathbf{P}=\mathbf{0}$, these nodes are ordered based on increasing values of associated 1-trees, adhering to the standard best-first rule.
In the neural version, when the disparity between the values of two subproblems is below a designated threshold, these subproblems are deemed equally meritorious. In such cases, the one exhibiting a superior expected value of the optimality score $O_T$ is prioritized.

§ GENERALIZATION TO ANY GRAPH DIMENSION

In the previous sections we assumed that the size of a GCN input matches the size of the TSP problem. The most straightforward approach to generalize this algorithm to accommodate graphs of any size is to train multiple models and select the appropriate one for the given instance size. However, this solution is practically infeasible due to the high computational cost of training a large number of networks with potentially small variations.

We propose here two distinct techniques that allow us to extract useful heuristic information to reach the exact optimal solution. This is achieved by using an approximation of the matrix $\mathbf{P}$ — at the expense, admittedly, of a degradation in the heuristic information conveyed in $\mathbf P$. Since we rely on the pre-trained GCN provided by <cit.> we use these techniques when the dimension of the input size is not $20$, $50$ or $100$, and so we do not have the neural network of the exact same size. Specifically, when the graph size is lower than the dimension of the utilized GCN, we suggest the "dummy cities" technique. Conversely, when the graph size exceeds $100$, the maximum dimension of the GCN, we suggest employing clustering. The following sub-sections detail these methods.

§.§ Dummy Cities

Let us assume to have a GCN model with input size $n$, when the size of the TSP problem is smaller than $n$, an effective strategy involves enlarging the instance graph with dummy cities until reaching the size $n$. These dummy cities are strategically placed near the original nodes, ensuring that the edges connecting them to the closest neighbor have a high probability of being included in the optimal solution. We use the augmented graph to obtain from the GCN the $p_{ij}$ estimations for all the known cities plus the dummy ones. Then we run <ref> on the original graph and using a $\textbf{P}$ restricted to the original cities. The intuition here is that the tour obtained from the optimal one on the augmented instance, after removing the dummy cities and connecting the remaining nodes, would closely approximate the real desired optimal tour. Similarly, we expect that the probability values $p_{ij}$ "obtained through" the GCN closely approximate those that would have been obtained if a model for such size was available. Empirical evidence for this intuition is provided in <ref>.

In this scenario, when we have the option to choose between enlarging or decreasing the number of cities, we expect that it is better to opt for the former. This is justified by the fact that the network can still analyze all the original edges and we can extract probabilities for each of them, instead of solely focusing on a subset of them.

§.§ Clustering

In case the TSP problem size is larger than the larger available GCN

we propose to preprocess the data in order to group nearby vertices into a number of clusters equal to the size that can be handled by the largest available GCN. Specifically, we utilize the  algorithm by  <cit.>  so that  vertices in the original problem are used as representatives of the clusters, and real edges can be then evaluated by the GCN. Once this is completed, we set the $p_{ij}$ equal to zero for of all the edges connecting cities within the same cluster and the other edges connecting two vertices from different clusters where at least one is not a medoid (i.e., we set to zero the weight of all edges that have not been explicitly seen by the neural network).

§.§ The Graph Convolutional Branch and Bound Algorithm

Our proposed algorithm GCBB, reported in <ref>, includes the steps necessary to handle problem instances where the number of cities does not match the size of one of the available networks. The algorithm essentially acts as a wrapper, initially leveraging the GCN including the dummy-city or the clustering preprocessing steps if they are necessary and then, upon obtaining the probability matrix $\mathbf{P}$, executing the BB <ref> to determine the optimal solution. More precisely for the pre-processing, if the condition at <ref> is met, the user must specify the location for adding the new dummy cities into the graph $G$ using the function at <ref>. Post-processing instead entails either removing rows and columns from $\mathbf{P}$ corresponding to the dummy cities or augmenting $\mathbf{P}$ with null probabilities for unseen edges.

Instance $G=(V,E)$, $\mathcal{M}$ the set of all the GCN models $M_d$ with different dimension $d$.
The optimal tour $T^\star$.

$\exists M_d \in  \mathcal{M}, \: d == |V|$
$\mathbf{P} = \texttt{graph\_conv\_net}(G, M_d)$

$d_{\max} = \max \{\:d\: \mid M_d \in \mathcal{M}\}$
$|V| < d_{\max} $
$d = \inf \{\:i\:\mid M_i \in \mathcal{M} \:\land \: i > |V|\} $
$\widetilde{G} = \texttt{add\_dummy\_cities}(G, d - |V| )$
$\widetilde{\mathbf{P}} = \texttt{graph\_conv\_net}(\widetilde{G}, M_{d})$
$\mathbf{P}= \texttt{reduce\_size}(\widetilde{\mathbf{P}})$

$\widetilde{\mathbf{P}} = \texttt{graph\_conv\_net}(\widetilde{G}, M_{d_{\max}})$
$P = \texttt{increase\_size}(\widetilde{\mathbf{P}})$

$T^\star = \texttt{branch\_and\_bound}(G,\mathbf{P})$*see <ref>
Graph Convolutional Branch and Bound

§ EXPERIMENTS

§.§ Graph instances

While the implemented branch and bound algorithm can be applied to various types of graph instances, the graph neural network is specifically trained on uniformly sampled Euclidean instances within the $[0,1]\times[0,1]$ square. Consequently, the experiments are ran through instances of this type whose cities are represented as pairs of coordinates $(x_i, y_i)$, and the travel distance between two nodes is the Euclidean distance. We note that every Euclidean graph can be rescaled to fit within the unit square without altering the set of optimal tours.

§.§ Experimental settings

Given the nature of floating-point numbers a threshold is defined below which two 1-trees are considered identical. To establish this threshold, a fraction of the upper bound proposed by <cit.> is selected. Ghosh proved that a valid upper bound for any random Euclidean instance in the unit square with $|V|=n$ vertices can be calculated as $1.27 \cdot \sqrt{n}$. In the implemented branch and bound we chose to consider two 1-trees as identical if they differ by less than $\texttt{EPSILON} = \frac{1.27 \cdot \sqrt{n}}{1000}$.

Moreover, given that we are tackling an NP-hard problem, it is reasonable to establish a time limit, beyond which, the computation is halted. When this happens the problem instance is deemed unsolved and its computation metrics disregarded. In our scenario, a  of 10 minutes is set, serving as the exit condition for the loop on <ref> of <ref>.

§.§ Comparison metrics

In our work we avoided the use of predefined libraries for standard functions because we need to monitor various metrics about the properties of the GCBB that otherwise would be hidden in the software. Before presenting the results obtained, we conclude this section by introducing the main metrics that we analyze later on.

In addition to conventional metrics like average solution time or average time to achieve optimality, we consider few metrics based on the number of nodes in the branch and bound tree, which has the advantage of being independent of the execution environment. The combinatorial explosion of these nodes is indeed at the core of the algorithm's high computational complexity. Specifically, we report:

* :  maximum depth of the generated tree;

* : depth at which the optimum is found;

* : total number of nodes required for completing the search for the optimum;

* : cardinality of the subset of Generated BB nodes necessitating a branching step;

* : number of nodes generated before reaching the optimum.

These metrics allow, for instance, to examine the contrast between  and ,  which gives an idea of the time saved by the PNN algorithm running on the heuristic information provided by the neural network.

We also report the value of the optimality score $\frac{\mathbb{E}[O_T^\star]}{n}$, which measures the goodness of the chosen criterion and so of the overall procedure. Indeed, we expect this value to be close to $1$ for all the 1-trees that represent the optimum of the analyzed problems.

§ RESULTS

We derive the value of the mentioned metrics by averaging it on the solutions obtained on $1000$ instances for all graph sizes.
In <ref>, we report results on instances where we had a model trained on the correct number of vertices. In <ref>, we report results where we had to resort to the techniques outlined in <ref> to accommodate for graph sizes for which a network of the correct size is not available.
In the tables, optimal values are highlighted in bold. If the disparities in the means are statistically significant according to a coupled Wilcoxon signed-rank test at the 5% significance level, the values are underlined.

The first three lines present the percentages of the total instances solved and how the optimum is found. In the next lines we report the mean values and their confidence intervals for all the chosen metrics. All these values are calculated when an exact GCN exists and refer to the same set of instances solved by both solvers within 10 minutes.

4c$n=20$     4c$n=50$     4c$n=100$

2cBB     2cGCBB     2cBB     2cGCBB     2cBB     2cGCBB

Instances solved     100.0%        100.0%        100.0%        100.0%        95.6%        95.6%

$\qquad$ with NN     8.9%        1.3%        0.0%        0.0%        0.0%        0.0%

$\qquad$ with PNN     -        86.0%        -       79.5%       -       48.12%

Total time (s)     0.01     0.0     1.55     0.0     1.65     0.85     3.25     0.77     51.69     5.91     45.73     5.4

BB time (s)     0.0     0.0      0.0     0.0     1.62     0.85     1.5     0.77     51.46     5.9     43.25     5.39

Time to Best (s)     0.0    0.0     0.0     0.0     0.75     0.36     0.51     0.24     33.6     4.09     19.1     3.21

BB tree depth      2.41      0.16      2.79      0.19     9.16     0.39      8.29      0.37     16.41      0.4      15.54     0.44

Depth of the optimum     1.61     0.08     0.36     0.07     5.83     0.21     1.65     0.22     11.03     0.26     6.44     0.43

Generated BB nodes      7.18      1.18     7.89     1.21      874.98     416.66     805.14     389.44      4439.6     500.27      3590.77     443.48

Explored BB nodes     6.34     1.05     7.46     1.02     778.19     398.87     753.88     376.0     3308.89     386.92     3013.26     380.32

BB nodes before optimum     3.58     0.65     1.61     0.68    388.68     191.23     214.31     127.41     2347.24     292.04     1561.14     273.65

$\frac{\mathbb{E}[O_{T^\star}]}{n}$    -       0.97     0.0    -          0.99     0.0    -       0.99     0.0

The first three lines present the percentages of the total instances solved and how the optimum is found. In the next lines we report the mean values and their confidence intervals for all the chosen metrics. All these values are calculated when an exact GCN does not exist and refer to the same set of instances solved by both solvers within 10 minutes.

4c$n=35$     4c$n=75$     4c$n=150$

2cBB     2cGCBB     2cBB     2cGCBB     2cBB     2cGCBB

Instances solved     100.0%        100.0%        99.2%        99.3%        43.6%        47.4%

$\qquad$ with NN     0.4%        0.2%        0.0%        0.0%        0.0%        0.0%

$\qquad$ with PNN     -        16.1%        -       4.43%       -       0.0%

Total time (s)     0.08     0.01     1.8     0.01     14.61     2.56     16.29     2.62     164.0     15.44     140.92     13.86

BB time (s)     0.06     0.01      0.06     0.01     14.51     2.55     13.88     2.61     163.37     15.4     137.94     13.83

Time to Best (s)     0.03    0.01     0.03     0.0     8.41     1.59     7.33     1.45     132.8     13.16     107.31     11.55

BB tree depth      5.86      0.3      5.82      0.31     13.57     0.41      13.43      0.41     16.93      0.41      17.25     0.43

Depth of the optimum     3.81     0.17     3.46     0.19     8.86     0.26     8.57     0.26     12.03     0.31     12.37     0.32

Generated BB nodes      88.02      14.51     83.72     13.91      2790.41     475.58     2619.97     474.98      2399.71     226.72      1995.95     202.76

Explored BB nodes     73.72     12.85     70.33     12.3     2264.99     402.53     2190.02     413.67     1483.21     143.73     1274.49     130.5

BB nodes before optimum     44.75     8.46     40.03     6.64    1455.35     278.99     1278.1     268.41     1527.24     166.92     1396.33     159.03

$\frac{\mathbb{E}[O_{T^\star}]}{n}$    -       0.67     0.01    -          0.85     0.0    -       0.38     0.0

As we can see in the two tables, GCBB achieves better results than classic BB when the graph size grows large. Indeed, the inference of the GCN requires polynomial time in the number of vertices, while the branch and bound is exponential, and, for large instances, the disparity between the two solvers widens as the neural network overhead gets amortized.

For graphs of $100$ vertices, the overall hybrid procedure requires less time on average for solving an instance, while for $n=20$ and $n=50$ the total time is smaller for the BB approach. It is worth noting however, that the metrics concerning branch and bound nodes were already significantly better with instances of size greater or equal than $35$ cities  (see <ref>). We expect this trend to emerge and even amplify with larger instances. Indeed, this is what we empirically found for instances with $150$ vertices, even if, in this scenario, we had to resort to non-optimal probabilistic values due to the lack of a GCN model able to solve instances with 150 vertices.

To better analyze the differences between the two solvers, we also display performance profiles, distinguishing between cases where an exact model exists (see <ref>) and cases where we had to adapt a model trained on instances of a different size (see <ref>). In each panel on the left, as the time threshold grows, we report the fraction of problem instances that are solved by each algorithm. Meanwhile, in the panel on the right, we consider also instances not yet solved. More precisely, as the allotted time grows, we report the proportion of instances for which each algorithm has reached a superior objective value. We note that, in these latter plots, all lines should approach zero if the algorithms are given enough time.

As we can observe, when the problem instance is small in size (e.g., 20 and 35 vertices), the overhead of the neural network (which requires few seconds to run) is not justified since all the problems are solved by BB in just few milliseconds. Instances with 50 and 75 cities are intermediate instances where the overhead is completely compensated by the quality of the extracted information, and a difference in times is not noticeable. With instances of 100 and 150 vertices, however, acquiring good heuristic information about the graph instance allows for a significant performance improvement. It's also worth noting that for graphs of these dimensions, not all instances are solved within the 10 minutes. In this situation, as it can be seen from the gap between the curves in <ref> for 100 and <ref> for 150 cities, the hybrid version terminates execution with a better solution in most cases.

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend style=at=(0.2,0.3),anchor=west,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/20_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/20_nodes/cumulative_profile.csv;

GCBB, BB;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,
legend style=at=(0.1,0.3),anchor=west,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/20_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/20_nodes/performance_profile.csv;

20 nodes.

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend pos=south east,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/50_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/50_nodes/cumulative_profile.csv;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/50_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/50_nodes/performance_profile.csv;

50 nodes.

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend pos=south east,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/100_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/100_nodes/cumulative_profile.csv;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/100_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/100_nodes/performance_profile.csv;

100 nodes.

Performances achieved by the GCBB and  BB solver grouped in two distinct performance profiles for each graph size where an exact model exists.

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend style=at=(0.2,0.3),anchor=west,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/35_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/35_nodes/cumulative_profile.csv;

GCBB, BB;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,
legend style=at=(0.1,0.3),anchor=west,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/35_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/35_nodes/performance_profile.csv;

35 nodes

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend pos=south east,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/75_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/75_nodes/cumulative_profile.csv;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/75_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/75_nodes/performance_profile.csv;

75 nodes

xlabel=Time (s),
ylabel=Proportion of instances solved,
major grid style = dotted,
legend pos=south east,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/150_nodes/cumulative_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/150_nodes/cumulative_profile.csv;

xlabel=Time (s),
ylabel=Proportion of instances with a better solution,
major grid style = dotted,

[no markers, color=nblue, line width=1.2pt] table [x=Time, y=Hybrid,col sep=comma]results/mean_results/150_nodes/performance_profile.csv;

[no markers, color=amber, line width=1.2pt] table [x=Time, y=Classic,col sep=comma]results/mean_results/150_nodes/performance_profile.csv;

150 nodes

Performances of the GCBB and BB solver grouped in distinct performance profiles for each graph size where an exact model does not exist. In panels <ref> and <ref>, the dummy cities technique is employed to increase the graph size, while clustering is used in <ref>.

xmin=0, xmax=150,
ymin=0, ymax=5000,
xtick = 20,35,50,75,100,150,
major grid style = dotted,
legend pos=north west,

[color=nblue,scatter, scatter src=explicit symbolic] coordinates (20, 7.890000)[a]

[color=amber,scatter, scatter src=explicit symbolic] coordinates (20, 7.184000)[a]

GCBB, BB

Let us conclude this section by analyzing the average value assumed by the optimality score $O_T$ on the optimal 1-trees. With the objective of maximizing this value across all 1-trees and the primary concern being differences in magnitude rather than exact values, we analyze the standardized score $\frac{\mathbb{E}[O_{T^\star}]}{n} \in [0,1]$, which provides a more intuitive interpretation of a score that assesses the heuristic proximity of the 1-tree to an optimal tour.
As we can see in <ref>, when we have a model whose size coincides with that of the instances we want to solve, the average score is close to 1, indicating that most of the edges predicted by the neural network as belonging to an optimal tour actually are.
In  <ref>, when we don't have a perfect model, the score understandably decreases, but it still remains significant, especially when we can add dummy cities, as mentioned in <ref>.
It is also worth noting that, in the case of clustering, this score significantly decreases since the neural network was only able to take into account $4,950$ edges out of the $11,175$ present in complete instances of 150 nodes (as we mentioned in <ref>, to avoid introducing bias, all probabilities of unobserved edges have been set to zero). Despite this, the hybrid solver outperforms the classical one in terms of the number of instances solved within 10 minutes. Also, it shows a decrease in runtime of approximately 25 seconds on average per instance.

§ CONCLUSION

Neural networks, as we currently know them, can not outperform traditional optimization approaches developed over the past decades for the exact resolution of NP problems. However, as this paper shows, it is possible to exploit the power of some deep learning systems to create hybrid algorithms to achieve better results. Indeed, there are many blind spots in classical optimization algorithms, where there is not a mathematically proven best thing to do among seemingly equivalent choices. To date, in these scenarios everyone has his own technique which he relies on because its validity has been empirically proven over the years. This is where the neural networks could help, as they can provide a good heuristic that can be trained using knowledge about the specific task at hand.
Potential future research directions for the traveling salesperson problem could involve the enhancement of state-of-the-art solvers, while, a more comprehensive approach, might involve the development of hybrid solvers to address other NP problems.

§ ACKNOWLEDGEMENT

This work was supported in part by the Spoke "FutureHPC & BigData" of the ICSC – Centro Nazionale di Ricerca in "High Performance Computing, Big Data and Quantum Computing", funded by European Union - NextGenerationEU; and in part by European Union within the H2020 RIA "European Processor Initiative — Specific Grant Agreement 2"[<https://www.european-processor-initiative.eu/>] under Grant 826647.