Finding Hamiltonian cycles with graph neural networks
We train a small message-passing graph neural network to predict Hamiltonian
cycles on Erd\H{o}s-R\'enyi random graphs in a critical regime. It outperforms
existing hand-crafted heuristics after about 2.5 hours of training on a single
GPU. Our findings encourage an alternative approach to solving computationally
demanding (NP-hard) problems arising in practice. Instead of devising a
heuristic by hand, one can train it end-to-end using a neural network. This has
several advantages. Firstly, it is relatively quick and requires little
problem-specific knowledge. Secondly, the network can adjust to the
distribution of training samples, improving the performance on the most
relevant problem instances. The model is trained using supervised learning on
artificially created problem instances; this training procedure does not use an
existing solver to produce the supervised signal. Finally, the model
generalizes well to larger graph sizes and retains reasonable performance even
on graphs eight times the original size.
Filip Bosnić, Mile Šikić
# Finding Hamiltonian cycles with graph neural networks

###### Abstract

We train a small message-passing graph neural network to predict Hamiltonian cycles on Erdos-Renyi random graphs in a critical regime. It outperforms existing hand-crafted heuristics after about 2.5 hours of training on a single GPU. Our findings encourage an alternative approach to solving computationally demanding (NP-hard) problems arising in practice. Instead of devising a heuristic by hand, one can train it end-to-end using a neural network. This has several advantages. Firstly, it is relatively quick and requires little problem-specific knowledge. Secondly, the network can adjust to the distribution of training samples, improving the performance on the most relevant problem instances. The model is trained using supervised learning on artificially created problem instances; this training procedure does not use an existing solver to produce the supervised signal. Finally, the model generalizes well to larger graph sizes and retains reasonable performance even on graphs eight times the original size.

 Machine learning, Neural nets, Graph algorithms, Heuristics design

## I Introduction

When dealing with problems that are computationally too costly to solve explicitly, such as NP-hard problems, it is common to rely on heuristics. The idea of using neural networks to train such heuristics is quite appealing and has attracted considerable interest over the years. One aims to enhance an algorithm, such as greedy search, with a neural network module that is trained to improve the decision-making of the algorithm. See [4, 8] or [29] for an introduction and an overview of the area. In practice, problem instances typically come from a distribution with specific biases which are hard to describe explicitly. These can be exploited by a neural network. As an illustration, let us consider the Hamiltonian cycle problem (HCP), which is at the core of this paper (nodes in the _cycle_ can not repeat). It asks the following:

**Problem 1** (HCP).: _Determine whether or not there exists a cycle that passes through all vertices of a given graph. If it exists, such a cycle is called a Hamiltonian cycle, and the graph is said to be Hamiltonian._

The general HCP is known to be NP-complete and thus computationally intractable. Currently, the fastest known exact solution algorithm is due to [5] and has worst-case complexity of \(\mathcal{O}(1.657^{n})\).

As far as applications are concerned, HCP is used to improve runtimes of rendering engines, see [2]. To do so, one solves the HCP for a dual graph of triangulation and renders the triangles in that order which reduces the number of points to process. Another application of HCP comes from genomics, more specifically, the problem of de novo genome assembly. The task here is to reconstruct the genetic material of an organism, i.e. the exact sequence of nucleobases on all of its chromosomes, from a large number of sequenced fragments called _reads_. As chromosomes contain hundreds of millions bases, correctly assembling a single one is already a huge undertaking, see [19] for an example. Interpreting overlaps between reads as edges, after preprocessing and cleaning (see [32]), one ends up with a _string graph_ as proposed in [20]. The Hamiltonian cycle in the string graph corresponds to the correct assembly of the chromosome. For more details see [22, 3, 28] and [14]. Both triangular meshes of 3d objects and string graphs of various assemblers (such as [3] or [28]) have specific structures and statistical properties arising from the context. These could make solving the HCP easier but are difficult to exploit directly. We show here how to exploit them using graph neural networks in a similarly specific setting of Erdos-Renyi random graphs.

For HCP in general, heuristics based on Hopfield networks were already trained in the early 90-ties, see [17, 18]. More recently, however, the area of geometric deep learning and graph neural networks has seen rapid developments and produced neural network layers such as message passing [9] or graph attention layers [30]. These layers are built to exploit any graph structure in data and can handle arbitrarily large graphs with a limited set of parameters, resembling convolution layers in computer vision. They have found applications in image and text processing, combinatorial optimization, physics, chemistry [9] and biology [22]. See [35] and [7] for a deeper dive into the area. In particular, they are excellent candidates for heuristics of graph-based problem. However, most efforts so far have been directed towards combinatorial optimization problems, the two-dimensional traveling salesman problem in particular. Heuristics for the 2d-TSP based on transformer architecture were trained in [16, 6] and those based on graph
neural networks in [34] and [12]. The state-of-the-art result is achieved in [6] where a comprehensive list of references can be found as well. It has to be noted that previously mentioned models still perform worse than the Concorde TSP solver [1], a state-of-the-art _exact_ solver based on branch and bound search combined with the cutting plane method. Nevertheless, theoretical complexities of neural network models are superior to Concorde. Let us also mention [13, 26] and [27] which work with general combinatorial optimization and constraint satisfaction problems.

In this paper we present a HCP solver based on _graph_ neural networks and show that it easily outperforms most hand-made heuristics. The code is available at [https://github.com/lbcb-sci/GNNs-Hamiltonian-cycles](https://github.com/lbcb-sci/GNNs-Hamiltonian-cycles).

## II Relation to TSP and 2d-TSP

It is known that the HCP can be reformulated as a special case of the _general traveling salesman problem (TSP)_:

**Problem 2** (TSP).: _Given a graph with a non negative length assigned to each edge, find the shortest cycle passing through all its vertices._

Hence, TSP solvers can be used for HCP and we shall exploit this by using _Concorde TSP solver_, see [1], to evaluate the performance of our models in Section V. While it is tempting to assume that all papers studying TSP are immediately applicable to the HCP, this _is not the case at all_. In particular, papers presenting neural network TSP solvers, such as [6, 12, 16] or [34] only study the special case of _two-dimensional TSP_:

**Problem 3** (2d-TSP).: _Given a set of points in the unit square \([0,1]^{2}\), find the shortest (in terms of Euclidean distance) cycle which passes through all of them._

The 2d-TSP introduces two simplifications to the general TSP:

* graphs are always fully connected and
* distances between nodes comply with Euclidean structure (triangle inequality).

Only \(2n\) point coordinates are required to describe a 2d-TSP instance, in contrast to \(n^{2}-n\) adjacency matrix weights needed for the general TSP. Moreover, 2d-TSP solvers can not be used to solve the HCP. On the contrary, we find it better to view the HCP and the 2d-TSP as two quite different aspects of the general TSP. The HCP focuses on complexities arising from discrete connectivity structure while the 2d-TSP deals with difficulties coming from the choice of edge lengths.

## III Problem setup

We only consider simple, undirected graphs and denote a typical graph example by \(G\) and its size (number of nodes) by \(n\). The HCP is classically posed as a decision problem: _Determine whether the graph contains a Hamiltonian cycle or not_. However, to put more emphasis on finding the actual cycle, which is important in practice, we also require that solvers produce at least one Hamiltonian cycle. In case the output of a solver is not a valid Hamiltonian cycle, which is straightforward to check, we assume the solver predicted that no Hamiltonian cycle exists.

### _Inputs and outputs_

A solver receives as input a graph \(G\) and outputs a walk \(v_{1}v_{2}\ldots v_{k}\) on \(G\) proposing a Hamiltonian cycle. The walk is considered to be closed if \(v_{1}=v_{k}\) and thus is a Hamiltonian cycle only if \(k=n+1\) and nodes \(v_{1},v_{2},\ldots v_{k-1}\) are all distinct.

### _Evaluation distribution_

The performance of HCP heuristics depends heavily on properties of graphs they are required to solve. Indeed, it is reasonable to have heuristics constructed specifically to achieve good performance on particular types of graphs, such as duals of triangulations or string graphs mentioned in Section I. As there are many possible applications of the HCP, finding a good class of evaluation graphs is a challenging task. Currently at least, there seems to be no agreed-upon class for this purpose. There are datasets of collected HCP problems, see, for example, [23] or [10], but they are not quite large enough to train neural networks on. A natural approach, used in early works such as [17, 18, 33] is to use random graphs generated by adding edges between pairs of vertices independently with a fixed probability \(p\in(0,1)\). Such random graphs are known as _Erdos-Renyi random graphs_ with edge probability \(p\). Papers working with 2d-TSP typically use a similar idea of evaluation on randomly generated problems, concretely the _random uniform euclidean (RUE)_ sets of two-dimensional points chosen uniformly at random from the unit square \([0,1]^{2}\).

However, using Erdos-Renyi graphs with _constant_ edge probability parameter \(p\) for evaluating the HCP has a major flaw. Intuitively it is clear that the HCP gets more difficult as the size of the graph increases. This is not the case for Erdos-Renyi graphs with _constant_\(p\) as indicated by Table I. It tracks performances of Concorde TSP solver and HybridHam heuristic from [25]. The performance of either solver clearly improves as the graph size increases, suggesting that the problem is in fact getting easier. The issue is that large graphs end up having too many edges, leading to many Hamiltonian cycles thus making it easier to find one.

This can be mended by carefully decreasing parameter \(p\) as the size of the graph increases. We rely on the following theorem from [15].

**Theorem 1** (Paraphrase of [15], Theorem 1.).: _Let \(\text{ER}(n,p)\) denote the Erdos-Renyi graph on \(n\) nodes with edge proba

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline  & \multicolumn{5}{c}{graph size} \\ \cline{2-6} Name & \(25\) & \(50\) & \(100\) & \(150\) & \(200\) \\ \hline Concorde & 0.80 & 1.0 & 1.0 & 1.0 & 1.0 \\ HybridHam & 0.41 & 0.68 & 0.79 & 0.84 & 0.87 \\ \hline \hline \end{tabular}
\end{table} TABLE I: Fraction of solved instances out of 5000 in supercritical regime, \(p=0.25\)
bility parameter \(p\). For every \(p_{H}\in(0,1)\) there is an explicit sequence \((p_{n})_{n\in\mathbb{N}}\) such that_

\[\mathbb{P}\left(\text{ER}(n,p_{n})\text{ is Hamiltonian}\right)\xrightarrow{n \rightarrow\infty}p_{H}.\]

_Concretely, one can take \(p_{n}=\frac{\ln n+\ln\ln n-\ln\ln p_{H}^{-1}}{n-1}\)._

In other words, for any \(p_{H}\) there is a procedure of generating graphs such that they contain a Hamiltonian cycle with a probability approximately equal to \(p_{H}\). We call this the _critical regime_ for the HCP. If the asymptotic behavior of \(p_{n}\) is above the one from the previous theorem, we speak _of the supercritical regime_. Examining the performance of Concorde solver in Table II shows that the empirical fraction of Hamiltonian cycles remains relatively stable and is fairly close to the asymptotic value of \(p_{H}=0.8\). By controlling the existence probability of Hamiltonian cycles we control their expected number in a graph and hence also the difficulty of the HCP. This motivates our use of Erdos-Renyi random graphs in the critical regime as the evaluation class. For simplicity, we use \(p_{H}=0.8\) for the rest of the paper although other values of \(p_{H}\) would work equally well. Two examples of random graphs in the critical regime are shown Fig. III.1.

### _Datasets_

We work exclusively with generated datasets. Our test dataset is sampled from the evaluation distribution described in the previous section and consists of \(5000\) Erdos-Renyi graphs in critical regime with \(p_{H}=0.8\) for each size \(n=25\), \(50\), \(100\), \(150\) and \(200\). This sample size is large enough so that the fraction of Hamiltonian graphs stays within \(\pm 2\%\) interval with \(95\%\) probability for every size \(n\). Train and validation datasets are generated from a different distribution described in Section IV-B. They are never explicitly sampled. Instead, graph examples are generated on the fly when needed. The train dataset is _limited_ to graphs of size \(25\) in order to emphasize generalization properties of the model.

## IV Model details

Our model is autoregressive and decodes the Hamiltonian cycle a single node at a time. It begins by selecting a starting node and then chooses between neighbors in each following step. The chosen node is then appended to the partial solution and the process repeats until a node gets visited twice. There are two main components, a _neural network component_ that guides the neighbor selection at each step and a _search algorithm_ which combines selected nodes into a Hamiltonian cycle. Concretely, given a _partial solution walk_\(v_{1}v_{2}\ldots v_{k}\) at \(k+1\)-th step of autoregressive decoding, the neural network component estimates with \(\mathcal{P}(v|v_{1}\ldots v_{k})\) the probability that extending the walk by node \(v\) will eventually lead to a Hamiltonian cycle (HC):

\[\mathcal{P}(v|v_{1}\ldots v_{k})\approx\mathbb{P}\left(v_{1}\ldots v_{k}v \subseteq\text{HC}\big{|}v_{1}\ldots v_{k}\subseteq\text{HC}\right).\]

The search algorithm then selects the neighbor \(v\) greedily according to estimated probabilities. It stops decoding when a node gets visited twice, i.e. \(v\in\{v_{1},\ldots v_{k}\}\), and returns \(v_{1}v_{2}\ldots v_{k}v\) as the solution. The greedy approach is the simplest case of beam search algorithm with beam width \(\beta=1\). For beam width \(\beta>1\), at each step \(k\) the algorithm keeps track of the top \(\beta\) partial walks according to score

\[\mathcal{S}(v_{1}v_{2}\ldots v_{k}) :=\prod_{j=1}^{k}\mathcal{P}(v_{j}|v_{1}\ldots v_{j-1})\] \[\approx\mathbb{P}(v_{1}v_{2}\ldots v_{k}\text{ is contained in a HC})\]

and extends them over all possible neighbors. A new set of top \(\beta\) partial solutions is then selected and the process repeats. Clearly, larger beam width \(\beta\) compensates for the performance of neural network at the cost of additional computations. While we report the performance of various beam widths in Table II, our basic model employs the simplest possible search algorithm (\(\beta=1\)) in order to emphasize the neural network part.

Our neural network uses _persistent node features_\(\mathbf{h}\) in the same way as in [31]. These features are passed on between applications of the neural network, adding a sort of recurrent structure to the network. This provides a way for the network to access information from previous decoding steps.

### _GNN architecture_

Since _graph neural networks (GNN)_ form the central component of our model, HCP information needs to be represented

Fig. 1: Examples of random graphs in the critical HCP regime. \(25\) nodes in top and \(50\) nodes in bottom row. Graphs in each row are identical. Right column graph is ordered in circle following a Concorde TSP solution, with the HP predicted by our basic model shown in solid red.
in the suitable form. We represent the adjacency matrix of the graph as a list of edges and one hot-encode the following three node-level feature channels. Two channels to mark the start and the end node of the partial solution plus a channel to mark all nodes the solution contains. Note that this is precisely the information needed to correctly extend the walk by an unvisited node or close the HC if necessary.

We employ the _encode-process-decode_ architecture analogue to the one used in [31]. This means that our GNN is divided into the _encoder_, _processor_ and _decoder_ networks. The whole GNN has around \(22\) thousand parameters. Both encoder and decoder are single layer, fully connected networks with ReLU activation that operate on node features _individually for each node_. The processor network, containing about \(95\%\) of all parameters, is the core part. It is a residual stack of \(5\) max-aggregation message passing layers, see [9] for more details. As names suggest, an input example is encoded, then processed and finally decoded by applying the above networks successively. In addition, we augment the output of the encoder with a randomized vector of features which was shown to improve the performance of GNNs in [24]. Algorithm 1 presents the pseudocode of a single forward pass. A "free" index \(i\in G\) in a line indicates that this line should be repeated for each node; symbol \(\bigoplus\) denotes concatenation in feature dimension; operator \(\max_{j\sim i}\) stands for maximum over neighbors of \(i\).

```
Input:\(G\) - graph with \(n\) vertices; \(\mathbf{x}\in\mathbb{R}^{(n,d_{n})}\) - partial walk repr.; \(\mathbf{h}\in\mathbb{R}^{(n,d_{n})}\) - persistent features Output:\(\mathbf{p}\in[0,1]^{n}\) - next-step probabilities per node. Hyperparams:\(d_{\text{in}}=3,d_{h}=32,d_{r}=4,n_{p}=5\) Params:\(\theta\equiv\{W_{E},b_{E},W_{P},b_{P},\ldots\}\) - NN weights // Encoder - Initialize features \(\mathbf{z}_{i}=W_{E}(\mathbf{x}_{i}\bigoplus\mathbf{h}_{i})+b_{E}\in\mathbb{ R}^{d_{h}-d_{r}}\) \(\mathbf{r}=\text{Uniform}\left([0,1]^{n\times d_{r}}\right)\in\mathbb{R}^{(n,d_ {r})}\) \(\mathbf{h}_{i}=\mathbf{z}_{i}\bigoplus\mathbf{r}_{i}\in\mathbb{R}^{d_{h}}\) // Processor - Apply residual max-MPNN layers for\(k=1,2,\ldots n_{p}\)do \(\mathbf{m}_{i}=\max_{j\sim i}\operatorname{ReLU}\left(W_{M}^{k}(\mathbf{h}_{i} \bigoplus\mathbf{h}_{j})+b_{M}\right)\in\mathbb{R}^{d_{h}}\) \(\mathbf{h}_{i}=\mathbf{h}_{i}+\operatorname{ReLU}\left(W_{P}^{(k)}\left( \mathbf{h}_{i}\bigoplus\mathbf{m}_{i}\right)+b_{P}^{(k)}\right)\in\mathbb{R}^ {d_{h}}\) // Decoder - Extract logits and probabilities \(\mathbf{l}_{i}=W_{D}(\mathbf{z}_{i}\bigoplus\mathbf{h}_{i})+b_{D}\in\mathbb{R}\) for\(i=1,2,\ldots,n\)do if\(i\sim\text{GetLastNode}(\mathbf{x})\)then \(\mathbf{l}_{i}=-\infty\) \(\mathbf{p}=\operatorname{softmax}\mathbf{l}\in\mathbb{R}^{n}\) return\(\mathbf{p}\), \(\mathbf{h}\)
```

**Algorithm 1**ApplyGNN\((G,\mathbf{x},\mathbf{h};\theta)\).

### _Training_

Our supervised approach requires a large number of solved HCP instances during training. Even though they can easily be generated using existing HCP solvers, we will show it is possible to train on artificially generated graphs such that HCP solution is known in advance. We believe that such methods are useful when working with problems similar to HCP for which no exact solvers are available. The construction of a training example starts from a graph \(G\) of arbitrary size but with no edges. A random permutation of nodes is then connected into a single cycle by adding appropriate edges into \(G\). This will be a Hamiltonian cycle in the final graph and is stored as a supervision signal. Finally, for every pair of vertices in \(G\) we add and edge connecting them with probability \(p_{\text{edge}}=0.125\) (independently of other pairs). \(p_{\text{edge}}\) is treated as a training hyperparameter and was determined through experimentation. While the distribution of training samples generated in this way is quite different from the evaluation distribution which consists of ER graphs, the results show that the basic model still generalizes well. Note also that the final graph may have Hamiltonian cycles other than the original one. All such cycles are ignored during training.

The training procedure is summarized in Algorithm 2. A single training example consists of a graph \(G\) and a Hamiltonian cycle \(v_{1}v_{2}\ldots v_{n}v_{1}\) on \(G\). The network is trained using _teacher forcing_ along this Hamiltonian cycle on the conditional cross-entropy loss \(\mathcal{L}\) defined by

\[\mathcal{L}\left(v_{1}\ldots v_{n}v_{1}\right)=-\sum_{i=2}^{n+1}\ln\left( \mathcal{P}(v_{i}|v_{1}\ldots v_{i-1})\right),\]

where \(v_{n+1}:=v_{1}\) for notational convenience. Remark that the summation index starts from \(2\) because the choice of the first node in a cycle is completely arbitrary. Loss \(\mathcal{L}\) is minimized over minibatches of 8 training examples using Adam optimizer with a learning rate of \(10^{-4}\) for 2000 epochs of 100 gradient updates. The final model checkpoint was selected based on the fraction of solved instances on validation set generated in the same way as the training set. The whole training was performed on a single NVIDIA GeForce RTX 3080 GPU and took about 2.5 hours. Weight initialization and other optimizer hyperparameters are kept to default PyTorch 1.11.0 values, [21].

## V Results and discussion

We evaluate the performance of our models by measuring the fraction of successfully solved problems on test dataset described in Section III and compared it with following heuristics:

1. _Concorde TSP solver_ - the state-of-the-art exact TSP solver from [1],
2. _HybridHam_ - an HCP heuristic from [25],
3. _Ant-inspired heuristic_ - an HCP heuristic presented in [33],
4. _Least degree first heuristic_ - simple greedy heuristic always selecting the neighbor with the lowest degree.

Let us remark that the ant-inspired heuristic is a convergence procedure which we terminate after \(5n^{2}\ln n\) steps. This bound matches the theoretical complexity of the basic model leading to a relatively fair comparison. In [33], authors suggest to
terminate after \(\mathcal{O}(n^{3})\) iterations but this is very time consuming. We list evaluation results in Table II and average inference times in Table III. Keeping in mind that testing can be performed on a different sample of \(5000\) graphs, the 95% confidence interval for all values in Table II is below \(\pm 0.02\). Models were run on a single NVIDIA GeForce RTX 3080 GPU while all other solvers were run on a single core of Intel Core i7-12700 processor. Note also that HybridHam, least degree first and ant-inspired heuristic were reimplemented in Python 3.8 and could be optimized for better performance.

Our HCP setup makes it impossible for a solver to produce a false positive prediction. Consequently, all solvers have perfect precision and metrics such as \(F_{1}\), \(F_{2}\) are unnecessarily complicated. As the number of true positives (solvable HCPs) is stable by construction of the evaluation set (0.8 in the limit), accuracy, recall and fraction of solved instances have similar qualitative behavior. Thus we only report the fraction of solved instances for each model.

In conclusion, after only a few hours of training our basic model clearly outperformed existing heuristic solvers without using any pre-solved HCP. We believe that techniques similar to the ones presented here can be used to quickly develop heuristic for variations or generalizations of the HCP. For example, the task of finding the longest cycle in a graph. Or the task of finding the route of minimal length which covers all the nodes in the graph (some of them maybe more than once). The class of Erdos-Renyi random graphs is used for simplicity and evaluation convenience since it allows for rough estimate of the difficulty of the HCP with respect to its size. Another class of graphs can be used just as well, provided that it is specific enough so that the neural network can exploit its statistical or structural peculiarities. But this typical happens with graph instances coming from practical problems. Moreover, polynomial complexity of \(\mathcal{O}(n^{2}\log n)\) for our basic model is superior to exponential complexity of exact solvers. For example, Concorde TSP solver on the RUE 2d-TSP instances was experimentally found to have complexity of \(\mathcal{O}(1.24\sqrt{n})\) in [11], although it is not clear how this translates to the critical regime HCP. Nevertheless, neural network solvers are yet to achieve reasonable performance on large input graphs and Concorde TSP solver remains the best-performing HCP solver. This comes as no surprise since Concorde also outperforms all existing neural network solvers for the 2d-TSP problem.

## VI Ablation study & training stability

The neural network component from Section IV is enhanced with persistent features and vectors of randomized features but can function without either of them. To estimate their importance, we separately removed each one and trained the corresponding reduced model 5 times from scratch. Average performances and confidence intervals of 2 standard deviation are shown in Fig. VI.1.

As shown on Fig. VI.1, persistent features play a crucial role in our model. Without them the model can fail to converge during training. This is probably because persistent features allow the model to updated its internal node representations throughout decoding process which results in an RNN-like behavior and consequently increases the range of message passing neural network layers. The use of randomized features is not as significant but becomes noticeable when generalizing to large graphs. Note also that Fig. VI.1 shows the standard deviation of training procedure for the main model to be around 5% of graphs solved.
