A multi-start local search algorithm for the Hamiltonian completion problem on undirected graphs
AA multi-start local search algorithm for the Hamiltonian completionproblem on undirected graphs
Jorik Jooken, Pieter Leyman, Patrick De Causmaecker
KU Leuven Kulak, Department of Computer Science, CODeS, Etienne Sabbelaan 53, 8500 Kortrijk (Belgium)[email protected], [email protected], [email protected]
Abstract
This paper proposes a local search algorithm for a specific combinatorial optimisation problem ingraph theory: the Hamiltonian Completion Problem (HCP) on undirected graphs. In this problem, theobjective is to add as few edges as possible to a given undirected graph in order to obtain a Hamilto-nian graph. This problem has mainly been studied in the context of various specific kinds of undirectedgraphs (e.g. trees, unicyclic graphs and series-parallel graphs). The proposed algorithm, however, con-centrates on solving HCP for general undirected graphs. It can be considered to belong to the categoryof matheuristics, because it integrates an exact linear time solution for trees into a local search algorithmfor general graphs. This integration makes use of the close relation between HCP and the minimum pathpartition problem, which makes the algorithm equally useful for solving the latter problem. Furthermore,a benchmark set of problem instances is constructed for demonstrating the quality of the proposed algo-rithm. A comparison with state-of-the-art solvers indicates that the proposed algorithm is able to achievehigh-quality results.
Keywords:
Metaheuristics; Matheuristics; Combinatorial optimisation; Hamiltonian completion prob-lem; Minimum path partition problem
The problem of determining whether a given undirected graph contains a Hamiltonian cycle is one of Karp’sfamous 21 NP-complete problems [21]. Being such a famous problem, it has been studied intensively formany years by the research community. This research has also led to many variations of this problem.One such variation is an interesting generalisation of the original problem in an optimisation context: theHamiltonian Completion Problem (HCP), which is the main topic of this paper.In this paper, we develop a matheuristic [8] for HCP on undirected graphs. We will identify a restrictedversion of this problem that can be solved exactly in polynomial time. The solution of this restricted problemis integrated into a multi-start local search algorithm that attempts to solve HCP on general undirectedgraphs. This algorithm will repeatedly perturb the solution of the restricted problem to obtain a high-qualitysolution for the original problem.The remainder of this paper is organised as follows: in section 2 a formal description of HCP is given.We further motivate this problem by discussing applications in section 3. Section 4 contains a brief overviewof earlier work that has been done for this problem. We highlight those results from literature that are relevantto understand the remainder of this paper. In section 5 the multi-start local search algorithm that we propose1 a r X i v : . [ c s . D M ] M a y or solving HCP is explained. The explanation is built in a bottom-up fashion: we first describe the separatecomponents of the algorithm, after which the complete algorithm is explained. In section 6 we will describehow a benchmark set of problem instances is generated. These problem instances will be used to comparethe algorithm from this paper with state-of-the-art algorithms. The results of these experiments are describedin section 7. Finally, a conclusion and possible avenues for further work are discussed in section 8. The Hamiltonian Completion Problem (HCP) on undirected graphs [17] is the following problem: given anundirected graph G = ( V , E ) , find a set of edges E (cid:48) such that G (cid:48) = ( V , E ∪ E (cid:48) ) is a Hamiltonian graph and thesize of E (cid:48) is as small as possible. A Hamiltonian graph is a graph that has a Hamiltonian cycle (a cycle thatvisits each node exactly once, except for the first and the last node, which are visited twice).As an example, consider Fig. 1. In this figure, the graph on the left is the graph G for which HCP hasto be solved. The graph in the middle shows G (cid:48) , which is obtained by adding 1 edge to G . The graph on theright illustrates that G (cid:48) is indeed a Hamiltonian graph (by showing the existence of a Hamiltonian cycle). Figure 1:
Graph G , graph G (cid:48) and a Hamiltonian cycle in G (cid:48) The need for solving HCP arises naturally in many applications. A first example application involvingfrequency assignment to transmitters has been identified in [13]. HCP can also be used to solve a special caseof the Travelling Salesman Problem (TSP) [4]. A problem instance of TSP on a complete edge-weightedgraph with only two different edge weights can be solved by solving HCP, as discussed in [26].Yet another application comes from a generalisation of the Bottleneck Travelling Salesman Problem[16]. In this generalisation, we are given an integer k and an edge-weighted undirected graph G = ( V , E ) .The goal is to find a set of at most k paths in G such that every node in V belongs to exactly one pathand such that the maximum edge weight amongst all edges in the paths is minimised. We assume that k ischosen such that at least one solution exists. This problem can be solved by doing a binary search for thesmallest maximum edge weight. To determine for a specific weight w whether a solution exists such thatthe maximum edge weight is at most w , one has to solve an instance of HCP. More specifically, let f ( G , w ) denote an unweighted, undirected graph obtained by removing from G those edges that have a larger weight2han w . The weights of the remaining edges are ignored to make the graph unweighted. The set of nodes in f ( G , w ) is the same as the set of nodes in G . Let E (cid:48) be the set of edges that is obtained by solving HCP on f ( G , w ) (i.e. adding the edges in E (cid:48) to f ( G , w ) results in a Hamiltonian graph). Now, there exists a solutionwith maximum edge weight at most w if and only if E (cid:48) contains at most k edges. Let w ∗ be the smallestvalue for which E (cid:48) contains at most k edges and let H be the graph obtained by adding the edges in E (cid:48) to f ( G , w ∗ ) (so H is Hamiltonian). Now the set of at most k paths in G (the answer to the problem) can befound by taking any Hamilonian cycle in H and dropping those edges that are in E (cid:48) . The required set ofpaths is simply the obtained set of connected components. All edges in the obtained paths will have a weightthat is at most equal to w ∗ . It is clear that HCP on undirected graphs is NP-hard, because it can be used to determine whether a graphhas a Hamiltonian cycle. Determining whether a graph has a Hamiltonian cycle, known as the Hamiltoniancycle problem, is NP-complete [15]. This means that HCP on undirected graphs cannot be solved exactly inpolynomial time, unless P = NP . Hence, most of the existing literature concentrates on efficiently solvingHCP for special cases of undirected graphs. Amongst those special cases are algorithms that solve HCP inpolynomial time for trees and unicyclic graphs [17][18][13], line graphs of trees [2][25], line graphs of cacti[10] and series-parallel graphs [27]. HCP has also been studied in the context of sparse random graphs [14].The amount of literature that concentrates on heuristically solving HCP is more limited [26][11][23].In the remainder of this section, we will introduce some definitions and results from existing literaturethat will be needed to understand the rest of the paper. Lemma 1 [11] states that HCP on undirected graphs can be solved by solving an instance of the symmetrictravelling salesman problem. Let G = ( V , E ) be an undirected, unweighted graph. Let G (cid:48) = ( V (cid:48) , E (cid:48) ) be anundirected, weighted, complete graph such that V = V (cid:48) and E (cid:48) consists of | V |∗ ( | V |− ) edges linking all pairsof nodes in V (cid:48) . The weight of an edge e ∈ E (cid:48) is defined to be equal to 0 if and only if there is also an edgein G that links the same pair of nodes and otherwise the weight is equal to 1. Lemma 1
HCP on G can be solved by determining the minimum weight Hamiltonian cycle in G (cid:48) andadding to G those edges of this cycle that have weight . Determining the minimum weight Hamiltoniancycle in G (cid:48) is an instance of the symmetric travelling salesman problem. Using this lemma, HCP can be solved by reducing it to an instance of the symmetric TSP problem andthen using state-of-the-art TSP solvers like Concorde [3] (an exact solver) or LKH [19] (a heuristic solver).When the graphs for which we want to solve HCP become bigger however, the computation time needed toproduce high-quality solutions for TSP with these solvers might no longer be acceptable.
In this subsection the minimum path partition problem is introduced. This problem is closely related to HCPand will be an important concept needed to understand the algorithm in section 5. First, some terminologyis introduced. 3 efinition 1
Two paths P and P are called vertex-disjoint if they do not have a vertex in common. Definition 2
A path partition of a graph G = ( V , E ) is a set of pairwise vertex-disjoint paths { P , P , . . . , P k } such that every node in V belongs to exactly path P i . Definition 3
A minimum path partition of graph G = ( V , E ) is a path partition of G, which consists of theleast amount of paths over all path partitions of G. Definition 4
The path partition number of a graph G (denoted PPN ( G ) ) is defined as the amount of pathsin a minimum path partition of G. Definition 5
The Hamiltonian completion number of a graph G (denoted HCN ( G ) ) is defined as the mini-mum amount of edges that need to be added to G to obtain a Hamiltonian graph. The minimum path partition problem on undirected graphs asks to find a minimum path partition on agiven undirected graph. The following lemma [17] shows that this problem is closely related to HCP onundirected graphs:
Lemma 2
Let G = ( V , E ) be an undirected graph. If G is Hamiltonian, then HCN ( G ) = and PPN ( G ) = .Else, HCN ( G ) = PPN ( G ) . In case
PPN ( G ) >
1, a minimum path partition of G can be used to determine which edges to add to G to make it Hamiltonian (by cyclically linking the endpoints of the paths of the path partition). This isillustrated in Fig. 2. The graph on the left of the figure is the graph G for which HCP needs to be solved.The minimum path partition is indicated on this graph. Nodes with an equal number are part of the samepath. The edges in green indicate the links between nodes that belong to the same path. Hence, the pathpartition number of the shown graph is equal to 3. The graph on the right of the figure illustrates that aHamiltonian graph can be obtained by cyclically linking the endpoints of the paths of the path partition. AHamiltonian cycle in this graph is indicated. Figure 2:
Graph G and a Hamiltonian cycle in graph G (cid:48) In case
PPN ( G ) =
1, either
HCN ( G ) = HCN ( G ) =
1, but there is no simple rule to deduct froma given minimum path partition which of these two equalities is the correct one. Note that
HCN ( G ) = PPN ( G ) = G such that there is4n edge between the endpoints of the single path in that minimum path partition. However, not necessarilyall minimum path partitions of G have this property. Hence, one cannot use the (non-)existence of an edgebetween the endpoints of a path to determine whether HCN ( G ) = HCN ( G ) = Lemma 3
Let G = ( V , E ) be an undirected graph that consists of multiple components (i.e. G is a discon-nected graph). HCP for G can be solved by computing a minimum path partition for every component of Gseparately. Using this lemma, it is clear that HCP for a graph G can be solved by solving HCP for every componentof G separately. Hence, in the rest of this paper we will assume that the graph G for which HCP has to besolved, is a connected graph, unless stated otherwise. If we restrict the input graph G , for which HCP has to be solved, to be a tree, we can efficiently solve HCP[18][13]: Lemma 4
Let G = ( V , E ) be a tree. It is possible to compute a minimum path partition of G with a timecomplexity of O ( | V | ) . The minimum path partition can then be used to solve HCP, as seen in the previous subsection. Thealgorithm to compute the minimum path partition of a tree is not further discussed in this paper, but interestedreaders are referred to [18] and [13].Furthermore, there is a relation between the PPN of a graph and the PPN of its spanning trees [17]:
Lemma 5
Let G = ( V , E ) be an undirected graph. For all spanning trees T of G, the following inequalityholds: PPN ( T ) ≥ PPN ( G ) . Furthermore there is at least one spanning tree T of G such that PPN ( T ) = PPN ( G ) . By combining lemma 4 and 5, it is clear that finding a minimum path partition for a graph can be doneby finding a minimum path partition for all of its spanning trees. Unfortunately, this approach is not alwaysfeasible, because the amount of spanning trees of a graph can be very large (e.g. a complete graph with n vertices has n n − spanning trees [9]). However, note that the fact that a graph has many spanning trees doesnot necessarily imply that it is hard to find a minimum path partition for this graph. As can be seen in theresults section, the algorithm that is proposed in this paper often is able to find the minimum path partitionfor graphs that have many spanning trees. The algorithm that we propose to solve HCP on an undirected graph G will attempt to find a minimum pathpartition of G . It does this indirectly, however, by trying to search for an optimal spanning tree T of G (i.e. such that PPN ( T ) = PPN ( G ) ). The algorithm starts with multiple initial spanning trees, that will beperturbed in order to obtain better spanning trees. A spanning tree T is called better than a spanning tree T if PPN ( T ) is smaller than or equal to PPN ( T ) . 5 .1 Initial spanning tree Several steps are taken to obtain an initial spanning tree of the graph G = ( V , E ) :Step 1: HCP is first transformed into a symmetric TSP instance, as discussed in section 4.1. A heuristicsolution for this symmetric TSP instance is then computed by the heuristic TSP solver LKH inorder to obtain an approximation for the minimum weight Hamiltonian cycle in the complete graph G (cid:48) (where G (cid:48) is the graph as described before lemma 1). In order to keep the running time low,the parameters that are supplied to LKH in this step are purposely chosen in such a way that onlya fraction of the full power of the LKH solver is being used (e.g. the supplied parameter RUNS ismuch smaller than the default amount of runs, see section 7.1).Step 2: This Hamiltonian cycle in G (cid:48) is used to obtain a path partition of G (by omitting those edges of thecycle that have weight 1).Step 3: This path partition is modified to obtain another path partition, by applying a rotation move [24] toall of its paths (in case it is possible to apply such a move). A rotation move modifies a path in thefollowing way: let P = ( v , v , . . . , v k ) be a path consisting of k vertices. If there is a vertex v i (with1 < i < k −
1) such that there is an edge between v i and v k in G , a rotation move can be appliedto obtain the path P (cid:48) = ( v , v , . . . , v i , v k , v k − , v k − , . . . , v i + ) . If there are multiple choices for v i , anarbitrary one is chosen. Fig. 3 illustrates the situation before step 3 is executed (on the left) andafterwards (on the right). A rotation move is applied for the first path (indicated with number 1). Figure 3:
The paths before and after applying rotation moves. The edges in green indicate the links between neigh-boring nodes in a path.
Step 4: The obtained path partition is transformed into a spanning tree of G by adding a set of edges of G , that do not already belong to the paths in the path partition. If the path partition consists of k paths, k − pre f erredRatiopre f erredRatio + (and from the category of normal edges withprobability 1 − pre f erredRatiopre f erredRatio + = pre f erredRatio + ). By defining the probabilities in this way, we havethe property that the parameter preferredRatio indicates the ratio between the probability that anedge is selected from the category of preferred edges and the probability that an edge is selectedfrom the category of normal edges. This parameter can be used to indicate how many times morelikely it is to select an edge from the category of preferred edges than an edge from of the categoryof normal edges. The spanning tree that is obtained, tends to have a lower path partition numberwhen more edges from the category of preferred edges are added. Hence, the intended use of thisparameter is to choose preferredRatio such that it is larger than 1, but this is not explicitly necessary.Algorithm 1 illustrates the pseudocode to obtain an initial spanning tree. Here, obtainInitialSpan-ningTree is a function that accepts as input a graph G and produces as output a spanning tree of G . Thefunction solveTSP is a function that accepts as input a graph and produces as output an approximation ofthe minimum weight Hamiltonian cycle in that graph. The function makeTree executes step 3 and step 4from above. It accepts as input a graph G and a path partition of G . It produces as output a spanning treeof G . The pseudocode of this function is given in algorithm 2. The function applyRotationMove accepts asinput a path p and returns as output the path that is obtained by applying a rotation move on p . Finally, thefunction addEdges corresponds to step 4 from above. It accepts as input a graph G and a path partition of G . It produces as output a spanning tree of G . Algorithm 1:
Obtain an initial spanning tree of the graph
Function obtainInitialSpanningTree( G ) Data:
G: an undirected graph
Result:
A spanning tree of G /* Start of code */
Let G (cid:48) = ( V (cid:48) , E (cid:48) ) be a complete graph such that V = V (cid:48) and the weight of an edge e ∈ E (cid:48) is definedto be equal to 0 if there is also an edge in G that links the same pair of nodes and 1 otherwise; /* step 1 */ minWeightHamiltonianCycle ← solveT SP ( G (cid:48) ) ; /* step 2 */ pathPartition ← /0; for edge in minWeightHamiltonianCycle :if weight ( edge ) = then pathPartition ← pathPartition ∪ edge ; /* step 3 and 4 */ return makeTree ( G , pathPartition ) To analyse the time complexity of obtaining an initial spanning tree, we will analyse the time complexityof the different steps separately. The time complexity of step 1 heavily depends on the parameters that are7 lgorithm 2:
Transform a path partition into a spanning tree
Function makeTree(
G, pathPartition ) /* Step 3 and 4 */
Data:
G: an undirected graphpathPartition: a path partition of G
Result:
A spanning tree of G /* Start of code */ newPathPartition ← pathPartition ; for path in newPathPartition :if rotation move can be applied on path then path ← applyRotationMove ( path ) ; return addEdges ( G , newPathPartition ) supplied to LKH. The parameters for LKH can be chosen such that the preprocessing step of this solver runsin nearly linear time (by choosing the POPMUSIC option for candidate set generation [20]) and the amountof local search steps is constant. For a more thorough treatment of the time complexity of this solver, thereader is referred to [19]. It is clear that step 2 and step 3 can each be implemented in O ( | V | + | E | ) . Inorder to efficiently execute step 4, one has to be able to quickly verify whether adding a certain edge wouldresult in a cycle (if this is the case, the edge should not be added). This is a well-known problem that canbe efficiently solved by using the union-find data structure [28]. By using this data structure, one can verifyin O ( α ( | V | )) time whether adding an edge would result in a cycle. Here, α ( | V | ) denotes the inverse ofthe Ackermann function [1]. Because of this, step 4 can be implemented in O ( | E | ∗ α ( | V | )) . Hence, thetotal time complexity to obtain an initial spanning tree is LKHTime ( V , E ) + O ( | V | + | E | ∗ α ( | V | )) , where LKHTime ( V , E ) denotes the time that the LKH solver spends in step 1. The spanning tree that is obtained, is repeatedly perturbed in order to obtain better spanning trees. A singleperturbation consists of multiple steps, which are explained below. In the following steps, we will denoteby G the graph for which HCP needs to be solved and by T a spanning tree of G that will be perturbed. Theprocess of perturbing a spanning tree can be thought of as marking some edges as active and other edges asinactive. Initially, all edges in G are inactive, except for the edges in T . After performing the perturbation,the set of active edges forms a (new) spanning tree of G . The steps of a single perturbation are explainedbelow.Step 1: Compute a minimum path partition of T . Mark all the edges that are part of T , but not part of thepath partition as inactive. Note that any path partition of T is also a path partition of G .Step 2: The main idea of this step is that we try to alter the path partition from the previous step to obtaina new path partition of G that consists of fewer paths. One possible way to approach this problemwould be to try to link together endpoints of the paths from the path partition by marking theconnecting edges as active. However, this idea can be generalised, without resulting in a highertime complexity: let PP = { P , P , . . . , P k } be the path partition obtained from step 1, which consistsof k paths. For all of the paths in PP , it is verified whether the path can be transformed into a cycle8y linking together the two endpoints of the path (i.e. it is verified whether there exists an edge in G such that the endpoints of this edge are equal to the endpoints of the path). In case the path canindeed be transformed into a cycle, this transformation is performed (the edge that is needed for thistransformation is marked as active). After doing this, we have a set of paths S (of size k ) and a setof cycles S (of size k ) such that k + k = k .These two sets will be iteratively updated by joining either a path with another path, a cycle with apath or a cycle with another cycle. To do this, the algorithm goes over all edges in G and verifieswhether the current edge can be used to join the structures (paths or cycles) at both endpoints of theedge. In case this is possible, the structures are joined and the sets S and S are updated. The threecases in which an edge can be used to join two structures are discussed next.Case 1: The edge is between two paths P i and P j . The edge can be used to join P i and P j if andonly if the edge is between an endpoint of P i and an endpoint of P j (and the paths P i and P j are not the same paths). If P i and P j can indeed be joined, the edge is marked as active,resulting in a longer path P (cid:48) . The set S is updated by removing P i and P j and adding P (cid:48) .Case 2: The edge is between a path P i and a cycle C j . The edge can be used to join P i and C j if andonly if exactly one endpoint of the edge corresponds to an endpoint of P i . If P i and C j canindeed be joined, the edge is marked as active. Furthermore, an edge of C j is marked asinactive such that the joined structure forms a path. There are two possible such edges thatcould be marked as inactive. An arbitrary one of these two edges is chosen. By joining P i and C j a new path P (cid:48) is obtained. The set S is updated by removing P i and adding P (cid:48) . Theset S is updated by removing C j .Fig. 4 illustrates case 2. The graph at the top of the figure consists of the path 1 − − − − − −
4. The green and blue edges represent the links betweenconsecutive nodes of the path and the cycle respectively. The black edge represents aninactive edge. There is an edge between the cycle and an endpoint of the path: the edgethat links node 3 and node 4. This edge is used to join the path and the cycle and willbe marked as active. There are two choices for an edge that could be marked as inactivein order to obtain a new path, namely the edge between node 4 and node 5 and the edgebetween node 4 and node 7. In this figure, the edge between node 4 and node 5 is markedas inactive in order to obtain the path 1 − − − − − −
5. The resulting path is shownat the bottom of the figure.Case 3: The edge is between a cycle C i and a cycle C j . The edge can be used to join the cycles ifand only if C i is not the same cycle as C j . In order to obtain a path by joining C i and C j ,the current edge is marked as active and furthermore one edge from C i and one edge from C j are marked as inactive. Again (as in case 2), there are two possible edges that can bemarked as inactive for a given cycle and an arbitrary one of these two is chosen for eachcycle. The set S is updated by adding the path P (cid:48) that results from joining the cycles. Theset S is updated by removing C i and C j .After having gone through all edges of G and joining the structures at the endpoints of the edgesin case this is possible, we have obtained a new set of paths S (cid:48) = { P , P , . . . , P k (cid:48) } and a new set ofcycles S (cid:48) = { C , C , . . . , C k (cid:48) } such that k (cid:48) + k (cid:48) ≤ k . Finally, for all cycles in S an arbitrary edgeis chosen to mark as inactive (to transform the cycle back into a path). The resulting set of activeedges forms a (new) path partition of G . 9 igure 4: A path and a cycle are joined to obtain a new path.
Step 3: This step is identical to step 3 from section 5.1. The path partition that is obtained from step 2 istransformed into another path partition by applying rotation moves to all of the paths of the pathpartition for which this is possible.Step 4: This step is identical to step 4 from section 5.1. The path partition from step 3 is transformed into aspanning tree T (cid:48) of G by marking a particular set of edges as active.To analyse the time complexity of perturbing a spanning tree, we will analyse the time complexity ofthe different steps separately. Recall from lemma 4 that step 1 can be implemented in O ( | V | ) . Step 2 can beimplemented in O ( | V | + | E | ) , because the state of an edge will only be changed at most once (from active toinactive or vice-versa). As we have seen in section 5.1, step 3 can be implemented in O ( | V | + | E | ) and step4 can be implemented in O ( | E | ∗ α ( | V | )) . Hence, the total time complexity of perturbing a spanning tree is O ( | V | + | E | ∗ α ( | V | )) . Now that the building blocks have been explained, we are ready to discuss the complete algorithm. Thealgorithm has three parameters: preferredRatio , maxAmountInitialSpanningTrees and maxAmountBadPer-turbations . The first parameter ( preferredRatio ) has already been discussed in step 4 of section 5.1, whereas maxAmountInitialSpanningTrees indicates how many initial spanning trees have to be generated. The thirdparameter ( maxAmountBadPerturbations ) requires a more detailled discussion of a crucial property of theperturbation operator. 10or every initial spanning tree T of G , a number of perturbations will be performed. Let perturbed ( T , G ) denote the spanning tree that is obtained by perturbing the spanning tree T (by making certain edges in G ac-tive and others inactive). The perturbation that was discussed before has the property that PPN ( perturbed ( T , G )) ≤ PPN ( T ) , so perturbing a spanning tree always yields another spanning tree that is at least as good. This prop-erty holds, because for every step of the perturbation, the path partition number of the graph that consists ofthe set of active edges does not increase. This property is the motivation for introducing the concept of a badperturbation of a spanning tree: a perturbation such that PPN ( perturbed ( T , G )) = PPN ( T ) . The parameter maxAmountBadPerturbations indicates how many bad perturbations are applied to a spanning tree, beforewe move on to the next spanning tree.For every spanning tree T (cid:48) , that is obtained after performing perturbations on an initial spanning tree T , a minimum path partition is computed. This path partition is used to estimate HCN ( G ) (see lemma 2).In case PPN ( T (cid:48) ) >
1, the algorithm estimates
HCN ( G ) to be equal to PPN ( T (cid:48) ) . In the special case wherethe path partition consists of a single path (i.e. PPN ( T (cid:48) ) = G between the endpoints of this path. The algorithm estimates HCN ( G ) as 1 in case such an edge does notexist and 0 otherwise.The pseudocode of the complete algorithm is shown in algorithm 3. Here, estimateHCN is a functionthat accepts as input a graph G for which HCP has to be solved and produces as output an estimate for HCN ( G ) . The function obtainInitialSpanningTree has been discussed in section 5.1. The function com-puteMinimumPathPartitionOfTree accepts as input a tree and produces as output a minimum path partitionof that tree. The function perturbed accepts as input a spanning tree of G as well as the graph G itself andproduces as output a (new) spanning tree of G (see section 5.2). Finally, the function estimateHCNPathPar-tition accepts as input a path partition of G as well as the graph G itself and produces as output an estimateof HCN ( G ) .To analyse the worst-case time complexity of the complete algorithm, first note that the path partitionnumber of any graph G = ( V , E ) is at most equal to | V | , because of the existence of a trivial path parti-tion in which every path consists of a single node. In every iteration of the while loop either the size orthe path partition decreases by at least one or the size of the path partition remains the same, in whichcase the amount of bad perturbations is incremented. If we combine this insight with the time complex-ities from the previous subsections, we obtain for the complete algorithm a worst-case time complexityof O ( maxAmountInitialSpanningTrees ) ∗ ( LKHTime ( V , E ) + ( | V | + maxAmountBadPerturbations ) ∗ ( | V | + | E | ∗ α ( | V | ))) . If the parameters are regarded as constants, the time complexity becomes LKHTime ( V , E ) + O ( | V | ∗ | E | ∗ α ( | V | )) . To the best of our knowledge, up until now there was no standard benchmark set of problem instancesalready available for HCP on undirected graphs. In order to be able to test the quality of the proposedalgorithm, we have constructed a benchmark set that consists of a variety of problem instances. A probleminstance for HCP on undirected graphs is simply an undirected graph for which HCP has to be solved.We have constructed a benchmark set of problem instances that vary with respect to size, density, degreedistribution and several other features. It is composed both of graphs that arise from practice as well asgraphs that were generated according to theoretical models. For the former case, a number of graphs wereselected from the instances of the tenth DIMACS challenge about graph partitioning and graph clustering[5]. Some of these graphs are directed graphs. These graphs are transformed into undirected graphs (byignoring the direction of the edges) to obtain the problem instances. The other problem instances were gen-11 lgorithm 3:
Estimate the HCN of a graph
Function estimateHCN( G ) Data:
G: an undirected graph
Result:
An estimate for
HCN ( G ) /* Start of code */ estimate ← ∞ ; for times in range( maxAmountInitialSpanningTrees ) : T ← obtainInitialSpanningTree ( G ) ; amountBadPerturbations ← pathPartition ← computeMinimumPathPartitionO f Tree ( T ) ; while amountBadPerturbations ≤ maxAmountBadPerturbations : T ← perturbed ( T , G ) ; newPathPartition ← computeMinimumPathPartitionO f Tree ( T ) ; if pathPartition and newPathPartition have the same amount of paths then amountBadPerturbations ← amountBadPerturbations + pathPartition ← newPathPartition ; newEstimate ← estimateHCNPathPartition ( pathPartition , G ) ; if newEstimate < estimate then estimate ← newEstimate ; if estimate = then /* Guaranteed global optimum */ return return estimate ;erated by using several graph generators. Some of these generators are part of the library of the StanfordNetwork Analysis Platform (SNAP) [22]. This library contains many graph generators that have an under-lying theoretical model. A brief overview of the six generators that were used to devise the benchmark setfollows.Generator 1: This generator follows the Erd˝os-R´enyi model [12]. This graph generator is a stochastic gen-erator and has two parameters n and p . The generated graph consists of n nodes and for everypair of nodes there is an edge between them with probability p . Hence, each (simple, undi-rected) graph that consists of n nodes and m edges is equally likely to be generated by thismodel (with probability p m ∗ ( − p ) n ∗ ( n − ) − m ).The problem instances that were generated by this generator contain a varying amount ofnodes (powers of 2 between 2 and 2 ). For a fixed amount of nodes n , several probleminstances were generated with varying probabilities (such that the expected average degree ofa node varies between 2 , , , . . . , n ).For graphs generated with this model, many theoretical results are known. We highlight tworelevant theoretical results from [7]: Theorem 1
Let < r < be a real number and let G ( n , p ) = ( V , E ) be a graph obtained by he Erd˝os-R´enyi random model with parameters n and p. If n tends to infinity and p is chosensuch that | E | = o ( | V | r ) , the probability that all components of the graph are trees, tends to 1. Because of the previous theorem, for sufficiently small values of p , the probability that thealgorithm proposed in this paper computes the global optimum tends to 1 if n tends to infinity.This is the case, because HCP is solved exactly for trees.The following theorem indicates for which graphs the Hamiltonian completion number islikely to be equal to 0: Theorem 2
Let f ( n ) be a function that tends to infinity if n tends to infinity and let G ( n , p ) =( V , E ) be a graph obtained by the Erd˝os-R´enyi random model with parameters n and p. Setp l = log ( n )+ log ( log ( n )) − f ( n ) n and p u = log ( n )+ log ( log ( n ))+ f ( n ) n . If n tends to infinity, the probabilitythat G ( n , p l ) contains a Hamiltonian cycle tends to 0 and the probability that G ( n , p u ) containsa Hamiltonian cycle tends to 1. Hence, for sufficiently large values of p the generated graphs are verly likely to have a Hamil-tonian completion number equal to 0.Generator 2: This generator generates circular graphs and has two parameters n and k . The generated graphconsists of n nodes, labelled 0 , , . . . , n −
1. For every node with label i , there is an edge to the k subsequent nodes ( i + ) mod n , ( i + ) mod n , . . . , ( i + k ) mod n . For all of these graphs, theHamiltonian completion number is equal to 0.The parameters that were used to generate problem instances by this generator are varying: n is between 500 and 30 ,
000 and k is between 3 and 10.Generator 3: This generator generates two-dimensional grid graphs and has two parameters n and m . Thegenerated graph consists of n ∗ m nodes, which can be thought of as n rows of m nodes whichare laid out in a two-dimensional grid. There is an edge between a node in row r and column c and a node in row r and column c if and only if | r − r | + | c − c | =
1. For these graphs,the Hamiltonian completion number is equal to 0 if either n or m are even and otherwise equalto 1.The problem instances that were generated by this generator contain between 300 and 10 , k other nodes is proportional with k − γ , where γ is a parameter ofthe generator (to be set by the user).The graphs that were generated according to this model contain between 300 and 3000 nodesand the nodes have an average degree varying between 6 and 16.Generator 5: This generator starts from a star graph consisting of n + n random edges toobtain a problem instance. A star graph consisting of n + K , n (a tree with 1 internal node and n leaves).The problem instances that were generated by this generator contain between 1000 and 4000nodes. 13enerator 6: The final generator generates structured trees and has two parameters l and c . The generatedgraphs are trees that consist of l levels of nodes and every node (except for the leave nodes) has c child nodes. Note that for these graphs, the algorithm proposed in the paper is guaranteed tocompute the global optimum, because our algorithm solves HCP exactly for trees.In order to not get results that are biased positively towards our own algorithm, only 2 probleminstances were generated by this generator (one instance was generated with l = c = l =
10 and c = In total, the benchmark set consists of 113 problem instances. This benchmark set is split into two parts:a tuning set (consisting of 50 randomly chosen problem instances) and a test set (consisting of the remain-ing 63 problem instances). Table 1 shows the distribution of the problem instances over each of the twosets. The tuning set will be used to find appropriate parameter settings for the algorithm proposed in thispaper. Furthermore, we will demonstrate the effect of varying these parameters and how this influencesthe performance of the algorithm. The test set will be used to show the added value of the algorithm com-ponents that were discussed in this paper. Furthermore, we compare our algorithm with the performanceof state-of-the-art algorithms and we attribute special attention to when our algorithm performs well andwhen it does not. All experiments were performed on a Intel i7-8550U CPU with a clock rate of 1.8GHz. The test data is available online at https://set.kuleuven.be/codes/hamiltoniancompletionproblem andat https://github.com/JorikJooken/ProblemInstancesHamiltonianCompletionProblem . Table 1:
Amount of problems in tuning setand test set
Tuning set Test setDIMACS 8 7Erd˝os-R´enyi 29 39circular 3 6grid 5 4preferential attachment 3 4star 2 1structured tree 0 2
As discussed before, the algorithm has three parameters: preferredRatio , maxAmountInitialSpanningTrees and maxAmountBadPerturbations . For the parameter preferredRatio , it is not immediately clear without anyexperiments which setting would lead to the best performance. For the other two parameters, however, it isclear that higher values correspond to a better performing algorithm, because they indicate the amount oflocal search that has to be done. Hence, for these two parameters we will try to find settings that yield agood trade-off between execution time and performance. Note that the LKH solver that is used in step 1 ofthe algorithm also has parameters itself. These parameters will not be subject to tuning, but instead they arechosen in such a way that only a fraction of the full power of the LKH solver is being used. The preciseparameters that were supplied to this solver can be found in table 2.14 able 2: Parameter settings for LKH for step 1 ofconstructing an initial spanning tree
Parameter Value
MOV E TY PE PATCHING C PATCHING A RUNS CANDIDAT E SET TY PE POPMUSICMAX T RIALS In the following sections, the multi-start local search algorithm from this paper will be used with variousparameter settings. For the sake of notational convenience, we will abbreviate the multi-start local searchalgorithm with a specific parameter setting as
MSLS a b c , where a , b and c correspond to the chosenvalues for preferredRatio , maxAmountInitialSpanningTrees and maxAmountBadPerturbations respectively.For every problem instance, the algorithm is interrupted after 1000 seconds if the algorithm has not finishedwithin this time. If the algorithm finishes in less than 1000 seconds, no further search is performed to fillthe remaining time. In order to obtain fair comparisons, in all of the experiments the reported averageswill be computed only over those problem instances for which all algorithms (in the scope of the specificexperiment) were able to produce an answer within 1000 seconds. Hence, the reported averages can varysignificantly between different experiments (depending on which problem instances the averages are com-puted over). This makes the reported averages only comparable within one specific experiment. We willreport separately about the problem instances for which an algorithm did not produce an answer within 1000seconds. All the execution times reported in this paper measure the time that passes between the start of thealgorithm and the end of the algorithm, unless stated otherwise. Note that this is different from the time thatpasses between the start of the algorithm and the first moment at which the reported result is found.In the first experiment, we have used several values (1, 5, 25 and 125) for the parameter preferredRatio .The summary of these results can be found in table 3. The first row of the table indicates on how many prob-lem instances the algorithm did not finish within 1000 seconds, whereas the second row of the table showsthe average amount of time needed to solve a problem instance. The third row displays the average amountof edges that the algorithm has added in order to obtain a Hamiltonian graph (the smaller the better). Fromthis table, we can see that varying this parameter has little effect (both regarding the amount of interruptions,the average execution time and the average amount of edges). The parameter setting of preferredRatio = Table 3:
Results for various values of preferredRatio
MSLS
MSLS
MSLS
25 10 3000
MSLS
125 10 3000amount interrupted 1 1 1 0average time (in seconds) 83 .
19 79 .
31 69 .
56 74 . .
53 234 .
51 234 .
51 234 . The second experiment that we have performed, concerns the effect of changing the parameter maxAm- untInitialSpanningTrees . The results of this experiment are summarised in table 4. From this table it isclear that increasing maxAmountInitialSpanningTrees does not improve the average amount of edges addeda lot. The average computing time does, however, increase quite fast. Note that the computing time ismultiplied by less than a factor of 10 if maxAmountInitialSpanningTrees is multiplied by 10. This happens,because the algorithm can stop considering new spanning trees as soon as a spanning tree with a HCN of 0 isfound. There are 8 problem instances for which the algorithm was not able to finish within 1000 seconds if maxAmountInitialSpanningTrees is equal to 100, but only 1 such problem instance if this parameter is equalto 10. This leads us to use a value of 10 for maxAmountInitialSpanningTrees . Table 4:
Results for various values of maxAmountInitialSpanningTrees
MSLS
25 1 3000
MSLS
25 10 3000
MSLS
25 100 3000amount interrupted 0 1 8average time (in seconds) 11 .
26 26 .
33 121 . .
10 39 .
00 38 . Finally, we performed an experiment where we have used several values (100, 300, 1000 and 3000) forthe parameter maxAmountBadPerturbations . The results of this experiment are summarised in table 5. Vary-ing this parameter has a larger impact on the result (average amount of edges added) than the previous twoparameters. By increasing this parameter, the amount of edges to add gets better (lower) and the executiontime gets worse (higher). Note, however, that again (as was the case for maxAmountInitialSpanningTrees )the ratio between the execution time and maxAmountBadPerturbations decreases if maxAmountBadPertur-bations increases. If we choose a value of 3000 for maxAmountBadPerturbations , there are quite someproblem instances for which the execution time is close to 1000 seconds. Hence, we have chosen to notincrease this parameter any further and we will use the value of 3000.
Table 5:
Results for various values of maxAmountBadPerturbations
MSLS
25 10 100
MSLS
25 10 300
MSLS
25 10 1000
MSLS
25 10 3000amount interrupted 0 1 1 1average time (in seconds) 11 .
52 20 .
89 44 .
48 69 . .
55 236 .
08 235 .
08 234 . These experiments lead us to use the following parameter configuration for testing purposes: preferre-dRatio = maxAmounInitalSpanningTrees =
10 and maxAmountBadPerturbations = The algorithm proposed in this paper makes use of the TSP solver LKH (in step 1 of the process of obtainingan initial spanning tree). In this subsection, we show that the other algorithm components contribute signif-icantly to the final result obtained by our algorithm. We make a comparison between the results obtainedby on the one hand
MSLS
25 10 3000 and on the other hand an algorithm that only performs step 1 of theprocess of obtaining an initial spanning tree (i.e.
MSLS
25 10 3000 without all other steps). We will call thelatter algorithm
MU LT I LKH . The same parameter settings were used for LKH in both algorithms. Theseparameter settings can be found in table 2. 16he results of running these two algorithms on the problem instances of the test set can be found in table6. From this table, we can see that about 20% of the total computing time of
MSLS
25 10 3000 was spenton step 1 for obtaining initial spanning trees. By also performing the other steps, the amount of added edgescould on average be reduced by approximately 7. Hence, the contribution of the other steps is significant.This is also confirmed by doing a paired samples t-test on the pairs obtained by pairing the results of bothalgorithms. The null hypothesis in this test states that the mean difference (the result of
MSLS
25 10 3000minus the result of
MU LT I LKH ) is equal to 0 and the alternative hypothesis states that the mean differenceis smaller than 0. A p-value of 0.004 was obtained, such that at a significance level of 1%, we reject the nullhypothesis in favor of the alternative hypothesis.
Table 6:
Results on test set for
MSLS
25 10 3000 and for
MULT I LKH
MSLS
25 10 3000
MULT I LKH amount interrupted 0 0average time (in seconds) 50 .
58 9 . .
44 165 . We will compare the results of our algorithm with the results obtained by the TSP solvers LKH and Concordeby regarding HCP as a specific TSP instance (see subsection 4.1). Concorde is an exact solver and hencefor this solver the most interesting aspect is its computation time. LKH on the other hand is a heuristicsolver. For LKH, we have used the recommended parameter settings [19] in this experiment. Note that theseparameter settings are much more time consuming than those used by LKH in step 1 of our own algorithm.The parameter settings used for LKH in this experiment can be found in table 7.
Table 7:
Recommended parameter settings for LKH
Parameter Value
MOV E TY PE PATCHING C PATCHING A RUNS CANDIDAT E SET TY PE POPMUSICMAX T RIALS | V | The results of running these three algorithms on the 63 problem instances of the test set are summarisedin table 8. All three algorithms were given a maximum computation time of 1000 seconds per test in-stance. The multi-start local search algorithm proposed in this paper needed on average around 72 .
7% and36 .
2% of the time needed by Concorde and LKH respectively. The average amount of edges added pro-duced by the three algorithms are very close to each other (Concorde produced the best results, followed by
MSLS
25 10 3000, followed by LKH). Recall that Concorde is an exact solver, so the result achieved byConcorde is the global optimum. From this comparison, we see that the other two algorithms are also oftenable to compute a solution that is (close to) the global optimum.17 able 8:
Results for Concorde,
MSLS
25 10 3000 and LKH
Concorde
MSLS
25 10 3000 LKHamount interrupted 4 0 6average time (in seconds) 42 .
38 30 .
82 85 . .
28 184 .
56 184 . There were 51 problem instances (out of 63) for which all three algorithms were able to produce theglobal optimum. The results on the other problem instances are shown in table 9. The values in the table arethe amount of edges added by each algorithm, whereas the corresponding computation times are providedbetween brackets.
MSLS
25 10 3000 was able to compute an answer within 1000 seconds for all probleminstances, while Concorde and LKH did not finish within this time limit for 4 and 6 problem instancesrespectively (indicated by the dashes in the table). The results obtained by
MSLS
25 10 3000 were at leastas good as those obtained by LKH for all problem instances except for one problem instance: a grid graphconsisting of 2 rows and 5000 columns ( grid graph preferential attachment preferential attachment
MSLS
25 10 3000 and LKH dif-fer significantly from the global optimum. The global optimum for preferential attachment
MSLS
25 10 3000 and LKH produced an answer of 9 and 12 respectively. Similarly the global optimum is 0for preferential attachment
MSLS
25 10 3000 and LKH produced an answer of 7 and 10 re-spectively. We suspect that the reason why these particular instances are hard to solve for
MSLS
25 10 3000and LKH is related to the degree distribution of the nodes. If the nodes in a graph all have a very low degree,the graph can be thought of as being quite similar to a tree, because the average degree of a node in a treeis approximately equal to 2. Hence, one could expect that for these graphs HCP can be solved (nearly) opti-mal, because HCP can be solved exactly in linear time for trees. If, however, the nodes in a graph all have avery high degree, the graph can be expected to contain many Hamiltonian cycles, such that the Hamiltoniancompletion number of these graphs will most likely be close to 0. This implies that one could also expectthat the HCP can be solved to (near-)optimality for these graphs. The degree distribution for the nodes inpreferential attachment graphs does, however, not belong to either of these two categories. Preferential at-tachment graphs have a power-law degree distribution and hence most nodes have a small degree while onlya few nodes have a very large degree in comparison with the other nodes.Finally, we feel the need to point out that it is possible to change the experimental setup such that thereexist parameter settings for the LKH solver that are able to obtain better results than the recommendedsettings. These improved parameter settings can be found in table 10. Recall that the execution timesmeasured in the rest of the paper measure the time that passes between the start of the algorithm and the endof the algorithm. By doing this, we give the algorithm the opportunity to be able to compute the best possibleanswer within the specified amount of time of 1000 seconds and hence the algorithm is not interrupted if noimprovement has been found after a long time. If we instead measure the time that passes between the startof the algorithm and the first moment upon which the reported answer is found, we obtain a different view.The results in this new experimental setting on the same problem instances as in table 9 can be found in table11. These results were kindly contributed to us by Keld Helsgaun, the creator of LKH. This experiment wasperformed on an iMac with a clock rate of 3.4 GHz. 18 able 9:
Results on problem instances for which at least one algorithm did not produce the global optimum.A dash indicates that the algorithm was not finished within 1000 seconds.
Name problem instance Concorde
MSLS
25 10 3000 LKH circle like − circle like − circle like − delaunay n − er
10 2 78 (6.31 s) 78 (40.81 s) 80 (26.9 s) er
13 11 0 (101.3 s) 0 (347.38 s) − er
13 3 12 (78.02 s) 12 (573.89 s) − grid graph − grid graph
50 50 − grid graph
80 80 − − preferential attachment preferential attachment Table 10:
Non-default parametersettings for LKH
Parameter Value
MAX CANDIDAT ES SUBGRADIENT NOREST RICT ED SEARCH NO
Table 11:
Results produced by LKH with non-defaultparameter settings. The execution times mea-sure the time that passes between the startof the algorithm and the first moment uponwhich the reported result is found. This is dif-ferent from the other times reported in this pa-per.
Name problem instance LKH circle like circle like circle like delaunay n
12 0 (2.04 s) er
10 2 78 (0.13 s) er
13 11 0 (8.00 s) er
13 3 12 (9.35 s) grid graph grid graph
50 50 0 (3.50 s) grid graph
80 80 0 (84.21 s) preferential attachment preferential attachment Conclusions and further work
In this paper, we have proposed a multi-start local search algorithm for HCP on undirected graphs. Further-more, we have constructed a benchmark set of problem instances for this problem, which can hopefully helpfuture researchers to be able to compare the performance of their algorithms against other algorithms. Thisbenchmark set was used to demonstrate that the algorithm in the paper performs well in comparison withstate-of-the-art TSP solvers, both in respect of time and quality of the produced answer.An interesting idea that was not further explored in this paper is related to decomposing a probleminstance into smaller problem instances. For HCP, the articulation points and bridges of the graph allow itto be decomposed, as discussed in [26]. It is, however, not immediately clear whether the overhead involvedin decomposing the problem and reassembling the different parts of the solution is worthwhile.An other interesting future research avenue is inspired upon the work that has been done for TSP. TheLKH TSP solver first converts the input graph to a new graph in which some edges, that are unlikely to bepart of an optimal TSP tour, are removed. By doing this, the search space can be drastically reduced and thesolver can focus its search on more promising parts of the search space. The most common measure that isused to estimate how likely it is that a certain edge is part of an optimal TSP tour is called the α -nearness ofan edge [19]. Note that for HCP the α -nearness of all edges are equal and hence this measure is not reallyuseful for HCP. It would be interesting to see whether some similar measure could be used for HCP. Acknowledgements
We gratefully acknowledge the support provided by the ORDinL project (FWO-SBO S007318N, DataDriven Logistics, 1/1/2018 - 31/12/2021). Pieter Leyman is a Postdoctoral Fellow of the Research Founda-tion - Flanders (FWO) with contract number 12P9419N.We also gratefully acknowledge the contribution of the data from table 10 and table 11 by Keld Helsgaun.
References [1] W. Ackermann. Zum hilbertschen aufbau der reellen zahlen.
Mathematische Annalen , 99(1):118–133, Dec 1928.[2] A. Agnetis, P. Detti, C. Meloni, and D. Pacciarelli. A linear algorithm for the hamiltonian completion number ofthe line graph of a tree.
Information Processing Letters , 79(1):17 – 24, 2001.[3] D. Applegate. Concorde - a code for solving traveling salesman problems. , 2001.[4] D. L. Applegate, R. E. Bixby, V. Chvatal, and W. J. Cook.
The Traveling Salesman Problem: A ComputationalStudy (Princeton Series in Applied Mathematics) . Princeton University Press, Princeton, NJ, USA, 2007.[5] D. A. Bader, H. Meyerhenke, P. Sanders, and D. Wagner, editors.
Graph Partitioning and Graph Clustering, 10thDIMACS Implementation Challenge Workshop, Georgia Institute of Technology, Atlanta, GA, USA, February 13-14, 2012. Proceedings , volume 588 of
Contemporary Mathematics . American Mathematical Society, 2013.[6] A.-L. Barab´asi and R. Albert. Emergence of scaling in random networks.
Science , 286(5439):509–512, 1999.[7] B. Bollob´as.
Modern Graph Theory , volume 184. 01 1998.[8] M. A. Boschetti, V. Maniezzo, M. Roffilli, and A. Boluf´e R¨ohler. Matheuristics: Optimization, simulationand control. In M. J. Blesa, C. Blum, L. Di Gaspero, A. Roli, M. Sampels, and A. Schaerf, editors,
HybridMetaheuristics , pages 171–177, Berlin, Heidelberg, 2009. Springer Berlin Heidelberg.
9] A. Cayley.
The Collected Mathematical Papers , volume 10 of
Cambridge Library Collection - Mathematics .Cambridge University Press, 2009.[10] P. Detti and C. Meloni. A linear algorithm for the hamiltonian completion number of the line graph of a cactus.
Discrete Applied Mathematics , 136(2):197 – 215, 2004. The 1st Cologne-Twente Workshop on Graphs andCombinatorial Optimization.[11] P. Detti, C. Meloni, and M. Pranzo. Local search algorithms for finding the hamiltonian completion number ofline graphs.
Annals of Operations Research , 156:5–24, 09 2007.[12] P. Erd¨os and A. R´enyi. On random graphs, I.
Publicationes Mathematicae (Debrecen) , 6:290–297, 1959.[13] D. S. Franzblau and A. Raychaudhuri. Optimal hamiltonian completions and path covers for trees, and a reduc-tion to maximum flow.
The ANZIAM Journal , 44(2):193204, 2002.[14] D. Gamarnik and M. Sviridenko. Hamiltonian completions of sparse random graphs.
Discrete Appl. Math. ,152(1-3):139–158, Nov. 2005.[15] M. R. Garey and D. S. Johnson.
Computers and Intractability; A Guide to the Theory of NP-Completeness . W.H. Freeman & Co., New York, NY, USA, 1990.[16] P. C. Gilmore and R. E. Gomory. Sequencing a one state-variable machine: A solvable case of the travelingsalesman problem.
Operations Research , 12(5):655–679, 1964.[17] S. Goodman and S. Hedetniemi. On the hamiltonian completion problem. In R. A. Bari and F. Harary, editors,
Graphs and Combinatorics , pages 262–272, Berlin, Heidelberg, 1974. Springer Berlin Heidelberg.[18] S. E. Goodman, S. T. Hedetniemi, and P. J. Slater. Advances on the hamiltonian completion problem.
J. ACM ,22(3):352–360, July 1975.[19] K. Helsgaun. An effective implementation of the linkernighan traveling salesman heuristic.
European Journalof Operational Research , 126(1):106 – 130, 2000.[20] K. Helsgaun.
Using POPMUSIC for Candidate Set Generation in the Lin-Kernighan-Helsgaun TSP Solver .Roskilde Universitet, 7 2018.[21] R. M. Karp.
Reducibility among Combinatorial Problems , pages 85–103. Springer US, Boston, MA, 1972.[22] J. Leskovec and R. Sosiˇc. SNAP: A general-purpose network analysis and graph-mining library.
ACM Transac-tions on Intelligent Systems and Technology (TIST) , 8(1):1, 2016.[23] H. P. Maretic and A. Grbic. A heuristics approach to hamiltonian completion problem (hcp). In , pages 1607–1612, May 2015.[24] L. P´osa. Hamiltonian circuits in random graphs.
Discrete Mathematics , 14(4):359 – 364, 1976.[25] A. Raychaudhuri. The total interval number of a tree and the hamiltonian completion number of its line graph.
Information Processing Letters , 56(6):299 – 306, 1995.[26] V. J. Rayward-Smith. The hamiltonian completion problem and its solution.
Engineering Optimization ,10(4):309–322, 1987.[27] K. Takamizawa, T. Nishizeki, and N. Saito. Linear-time computability of combinatorial problems on series-parallel graphs.
J. ACM , 29(3):623–641, July 1982.[28] R. E. Tarjan and J. van Leeuwen. Worst-case analysis of set union algorithms.
J. ACM , 31(2):245–281, Mar.1984., 31(2):245–281, Mar.1984.