Generating subtour elimination constraints

for the TSP from pure integer solutions

Ulrich Pferschy∗ Rostislav Staněk∗
arXiv:1511.03533v1 [math.OC] 11 Nov 2015

Abstract
The traveling salesman problem (TSP) is one of the most prominent combinat-
orial optimization problems. Given a complete graph G = (V, E) and non-negative
distances d for every edge, the TSP asks for a shortest tour through all vertices with
respect to the distances d. The method of choice for solving the TSP to optimality
is a branch and cut approach. Usually the integrality constraints are relaxed first
and all separation processes to identify violated inequalities are done on fractional
solutions.
In our approach we try to exploit the impressive performance of current ILP-
solvers and work only with integer solutions without ever interfering with fractional
solutions. We stick to a very simple ILP-model and relax the subtour elimination
constraints only. The resulting problem is solved to integer optimality, violated
constraints (which are trivial to find) are added and the process is repeated until a
feasible solution is found.
In order to speed up the algorithm we pursue several attempts to find as many
relevant subtours as possible. These attempts are based on the clustering of ver-
tices with additional insights gained from empirical observations and random graph
theory. Computational results are performed on test instances taken from the
TSPLIB95 and on random Euclidean graphs.

Keywords. traveling salesman problem; subtour elimination constraint; ILP solver; ran-
dom Euclidean graph

1 Introduction
The Traveling Salesman/Salesperson Problem TSP is one of the best known and most
widely investigated combinatorial optimization problems with four famous books entirely
devoted to its study ([11], [18], [8], [2]). Thus, we will refrain from giving extensive
references but mainly refer to the treatment in [2]. Given a complete graph G = (V, E)
with |V | = n and |E| = m = n(n − 1)/2, and nonnegative distances de for each e ∈ E,
the TSP asks for a shortest tour with respect to the distances de containing each vertex
exactly once.
∗
{pferschy, rostislav.stanek}@uni-graz.at. Department of Statistics and Operations Research,
University of Graz, Universitätsstraße 15, 8010 Graz, Austria

1
Let δ(v) := {e = (v, u) ∈ E | u ∈ V } denote the set of all edges adjacent to v ∈ V .
Introducing binary variables xe for the possible inclusion of any edge e ∈ E in the tour
we get the following classical ILP formulation:
X
minimize de xe (1)
e∈E
X
s.t. xe = 2 ∀ v ∈ V, (2)
e∈δ(v)
X
xe ≤ |S| − 1 ∀ S ⊂ V, S 6= ∅, (3)
e=(u,v)∈E
u,v∈S

xe ∈ {0, 1} ∀e∈E (4)

(1) defines the objective function, (2) is the degree equation for each vertex, (3) are the
subtour elimination constraints, which forbid solutions consisting of several disconnected
tours, and finally (4) defines the integrality constraints. Note also that some subtour
elimination constraints are redundant: For the vertex sets S ⊂ V , S 6= ∅, and S ′ = V \S
we get pairs of subtour elimination constraints both enforcing the connection of S and
S ′.
The established standard approach to solve TSP to optimality, as pursued successfully
during the last 30+ years, is a branch-and-cut approach, which solves the LP-relaxation
obtained by relaxing the integrality constraints (4) into xe ∈ [0, 1]. In each iteration of
the underlying branch-and-bound scheme cutting planes are generated, i.e. constraints
that are violated by the current fractional solution, but not necessarily by any feasible
integer solution. Since there exists an exponential number of subsets S ⊂ V implying
subtour elimination constraints (3), the computation starts with a small collection of
subsets S ⊂ V (or none at all), and identifies violated subtour elimination constraints
as cutting planes in the so-called separation problem. Moreover, a wide range of other
cutting plane families were developed in the literature together with heuristic and exact
algorithms to find them (see e.g. [20, ch. 58], [2]). Also the undisputed champion among
all TSP codes, the famous Concorde package (see [2]), is based on this principle.
In this paper we introduce and examine another concept for solving the TSP. In
Section 2 we introduce the basic idea of our approach. Some improvement strategies
follow in Section 3 with our best approach presented in Subsection 3.6. Since the main
contribution of this paper are computational experiments, we discuss them in detail
in Section 4. The common details of all these tests will be given in Subsection 4.1. In
Section 5, we present some theoretical results and further empirical observations. Finally,
we provide an Appendix with illustrations, graphs and two summarizing tables (Tables 8
and 9).

2 General solution approach

Clearly, the performance of the above branch-and-cut approach depends crucially on the
performance of the used LP-solver. Highly efficient LP-solvers have been available for

2
quite some time, but also ILP-solvers have improved rapidly during the last decades and
reached an impressive performance. This motivated the idea of a very simple approach
for solving TSP without using LP-relaxations explicitly.
The general approach works as follows (see Algorithm 1). First, we relax all subtour
elimination constraints (3) from the model and solve the remaining ILP model (cor-
responding to a weighted 2-matching problem). Then we check if the obtained integer
solution contains subtours. If not, the solution is an optimal TSP tour. Otherwise, we
find all subtours in the integral solution (which can be done by a simple scan) and add
the corresponding subtour elimination constraints to the model, each of them represen-
ted by the subset of vertices in the corresponding subtour. The resulting enlarged ILP
model is solved again to optimality. Iterating this process clearly leads to an optimal
TSP tour.
Input: TSP instance
Output: an optimal TSP tour
1: define current model as (1), (2), (4);
2: repeat
3: solve the current model to optimality by an ILP-solver;
4: if solution contains no subtour then
5: set the solution as optimal tour;
6: else
7: find all subtours of the solution and add the corresponding subtour elimination
constraints into the model;
8: end if
9: until optimal tour found;

Algorithm 1: Main idea of our approach.

Every execution of the ILP-solver (see line 3) will be called an iteration. We define
the set of violated subtour elimination constraints as the set of all included subtour
elimination constraints which were violated in an iteration (see line 7). Figures 5 and
8 – 19 in the Appendix illustrate a problem instance and the execution of the algorithm
on this instance respectively.
It should be pointed out that the main motivation of this framework is its simplicity.
The separation of subtour elimination constraints for fractional solutions amounts to the
solution of a max-flow or min-cut problem. Based on the procedure by Padberg and
Rinaldi [15], extensive work has been done to construct elaborated algorithms for per-
forming this task efficiently. On the contrary, violated subtour elimination constraints of
integer solutions are trivial to find. Moreover, we refrain from using any other additional
inequalities known for classical branch-and-cut algorithms, which might also be used to
speed up our approach, since we want to underline the strength of modern ILP-solvers
in connection with a refined subtour selection process (see Section 3.6).
This approach for solving TSP is clearly not new but was available since the earliest
ILP formulation going back to [6] and can be seen as folklore nowadays. Several authors

3
followed the concept of generating integer solutions for some kind of relaxation of an
ILP formulation and iteratively adding violated integer subtour elimination constraints.
However, it seems that the lack of fast ILP-solvers prohibited its direct application in
computational studies although it was used in an artistic context by [4].
Miliotis [12] also concentrated on generating integer subtour elimination constraints,
but within a fractional LP framework. The classical paper by Crowder and Padberg [5]
applies the iterative generation of integer subtour elimination constraints as a second part
of their algorithm after generating fractional cutting planes in the first part to strengthen
the LP-relaxation. They report that not more than three iterations of the ILP-solver
for the strengthened model were necessary for test instances up to 318 vertices. Also
Grötschel and Holland [7] follow this direction of first improving the LP-model as much
as possible, e.g. by running preprocessing, fixing certain variables and strengthening the
LP-relaxation by different families of cutting planes, before generating integer subtours
as last step to find an optimal tour. It turns out that about half of their test instances
never reach this last phase. In contrast, we stick to the pure ILP-formulation without
any previous modifications.
From a theoretical perspective, the generation of subtours involves a certain trade-off.
For an instance (G, d) there exists a minimal set of subtours S ∗ , such that the ILP model
with only those subtour elimination constraints implied by S ∗ yields an overall feasible,
and thus optimal solution. However, in practice we can only find collections of subtours
larger than S ∗ by adding subtours in every iteration until we reach optimality. Thus,
we can either collect as many subtours as possible in each iteration, which may decrease
the number of iterations but increases the running time of the ILP-solver because of the
larger number of constraints. Or we try to control the number of subtour elimination
constraints added to the model by trying to judge their relevance and possibly remove
some of them later, which keeps the ILP-model smaller but may increase the number of
iterations. In the following we describe various strategies to find the “right” subtours.

2.1 Representation of subtour elimination constraints

The subtour elimination constraints (3) can be expressed equivalently by the following
cut constraints:
X
xe ≥ 2 ∀ S ⊂ V, S 6= ∅ (5)
e=(u,v)∈E
u∈S,v6∈S

Although mathematically equivalent, the two ways of forbidding a subtour in S may

result in quite different performances of the ILP-solver.
It was observed that in general the running time for solving an ILP increases with
the number of non-zero entries of the constraint matrix. Hence, we also tested a hybrid
variant which chooses between (3) and (5) by picking for each considered set S the version

4
with the smaller number of nonnegative coefficients on the left-hand side as follows:

if |S| ≤ 2n+1
P
e=(u,v)∈E xe ≤ |S| − 1 3
∀ S ⊂ V, S 6
= ∅
P u,v∈S (6)
e=(u,v)∈E xe ≥ 2 if |S| > 2n+1
3
u∈S,v6∈S

We performed computational tests of our approach to compare the three representations

of subtour elimination constraints, namely (3), (5) and (6), and list the results in Table 1.
Technical details about the setup of the experiments can be found in Subsection 4.1.

instance s.e.c. as in (3) s.e.c. as in (5) s.e.c. as in (6)

sec. #i. #c. sec. #i. #c. sec. #i. #c.
kroA150 89 12 82 75 12 82 62 12 82
kroB150 52 13 77 237 13 77 54 13 77
u159 9 5 39 13 5 39 9 5 38
brg180 62 14 56 36 5 29 64 16 67
kroA200 2153 11 95 1833 11 95 2440 11 95
kroB200 45 7 65 146 7 65 37 7 65
tsp225 149 15 102 376 16 105 155 16 106
a280 114 10 59 249 10 56 132 10 63
lin318 7171 13 177 8201 13 177 7158 13 177
gr431 5973 22 186 19111 22 187 5925 22 186
pcb442 4406 43 215 6186 41 197 2393 43 207
gr666 33259 14 216 189421 14 217 40111 14 216
mean ratio (sec.) 2.305960 0.971694
RE A 150 23 12 61 65 12 61 26 12 61
RE A 200 81 15 84 139 15 84 76 15 84
RE A 250 156 14 82 208 14 82 133 14 82
RE A 300 534 14 123 4819 14 123 692 14 123
RE A 350 404 9 110 789 9 110 650 9 110
RE A 400 49234 16 179 247511 16 179 24619 16 179
RE A 450 4666 8 117 13806 8 117 3022 8 117
RE A 500 68215 12 167 155977 12 167 30809 12 167
mean ratio (sec.) 3.390678 0.928176
mean ratio all 2.739847 0.954287

Table 1: Comparison of the behavior of the algorithm for different representations of

subtour elimination constraints. Mean ratios refer to the arithmetic means over ratios
between the running times of the approaches using the subtour elimination constraints
represented as in (5) and (6) respectively and the running time of the approach using
the subtour elimination constraints represented as in (3). “sec.” is the time in seconds,
“#i.” the number of iterations and “#c.” the number of subtour elimination constraints
added to the ILP before starting the last iteration.

5
It turned out that the three versions sometimes (but not always) lead to huge dif-
ferences in running time (up to a factor of 5). This is an interesting experience that
should be taken into consideration also in other computational studies. From our lim-
ited experiments it could be seen that version (5) was inferior most of the times (with
sometimes huge deviations) whereas only a small dominance of the hybrid variant (6) in
comparison with the standard version (3) could be observed. This is due to the small size
of most subtours occurring during the solution process (the representation (3) equals to
the representation (6) in these cases). But since also bigger subtours can occur (mostly
in the last iterations), we use the representation (6) for all further computational tests.
For more details about different ILP-models see [14].

3 Generation of subtours
As pointed out above, the focus of our attention lies in the generation and selection of
a “good” set of subtour elimination constraints, including as many as possible of those
required by the ILP-solver to determine an optimal solution which is also feasible for
TSP, but as few as possible of all others which only slow down the performance of the
ILP-solver.
Trying to strike a balance between these two goals we followed several directions,
some of them motivated by theoretical results, others by visually studying plots of all
subtours generated during the execution of Algorithm 1.

3.1 Subtour elimination constraints from suboptimal integer solutions

Many ILP-solvers report all feasible integer solutions found during the underlying branch-
and-bound process. In this case, we can also add all corresponding subtour elimination
constraints to the model. These constraints can be considered simply as part of the set of
violated subtour elimination constraints. Not surprisingly, these additional constraints
always lead to a decrease in the number of iterations for the overall computation and to an
increase in the total number of subtour elimination constraints generated before reaching
optimality (see Table 2). While the time consumed in each iteration is likely to increase,
it can also be observed that the overall running time is often decreased significantly by
adding all detected subtours to the model. On the other hand, for the smaller number of
instances where this is not the case, only relatively modest increases of running times are
incurred. Therefore, we stick to adding all detected subtour elimination constraints for
the remainder of the paper. The algorithm in this form will be called BasicIntegerTSP.

3.2 Subtours of size 3

The next idea we tried was to add subtour inequalities corresponding to some subtours of
size 3 into the model before starting the iteration process (i.e. in line 1 of Algorithm 1).
This idea was motivated by the observations that in many examples smaller subtours
(with respect to their cardinality) occur more often than the larger ones. However, there
are |V3 | such subtours and thus we should concentrate only on a relevant subset of them.

6
only subtours all subtours:
from ILP-optima BasicIntegerTSP
instance sec. #i. #c. sec. #i. #c.
kroA150 62 12 82 19 7 136
kroB150 54 13 77 179 8 148
u159 9 5 38 6 4 49
brg180 64 16 67 44 4 103
kroA200 2440 11 95 677 8 237
kroB200 37 7 65 31 5 121
tsp225 155 16 106 178 9 261
a280 132 10 63 157 11 143
lin318 7158 13 177 6885 8 357
gr431 5925 22 186 2239 9 453
pcb442 2393 43 207 2737 11 501
gr666 40111 14 216 17711 8 789
mean ratio (sec.) 0.946130
RE A 150 26 12 61 23 8 100
RE A 200 76 15 84 72 7 163
RE A 250 133 14 82 138 9 186
RE A 300 692 14 123 866 6 295
RE A 350 650 9 110 411 5 252
RE A 400 24619 16 179 8456 8 454
RE A 450 3022 8 117 2107 5 279
RE A 500 30809 12 167 15330 6 436
mean ratio (sec.) 0.786451
mean ratio all 0.882259

Table 2: Using all constraints generated from all feasible solutions found during the
solving process vs. using only the constraints generated from the final ILP solutions of
each iteration. Mean ratios refer to the arithmetic means over ratios between the running
times of BasicIntegerTSP over the other approach. “sec.” is the time in seconds, “#i.”
the number of iterations and “#c.” the number of subtour elimination constraints added
to the ILP before starting the last iteration.

After studying our computational tests we decided to use the shortest ones with respect
to their length. Table 3 summarizes our computational results and it can be seen that
this idea actually tends to slow down our approach. Thus we did not follow it any more.

3.3 Subtour selections

As mentioned above, a large number of subtour inequalities which are not really needed
only slow down our approach. Thus we also tried not to use all subtour inequalities we

7
1 1
instance p=0 p = 10000 p = 1000
sec. #i. #c. sec. #i. #c. sec. #i. #c.
kroA150 19 7 136 19 7 97 40 5 116
kroB150 179 8 148 71 7 178 134 5 105
u159 6 4 49 8 4 46 6 3 24
brg180 44 4 103 34 15 108 82 9 270
kroA200 677 8 237 879 5 157 504 4 133
kroB200 31 5 121 32 5 61 43 5 60
tsp225 178 9 261 149 10 224 167 9 202
a280 157 11 143 138 9 98 156 6 101
lin318 6885 8 357 5360 8 302 1435 8 291
gr431 2239 9 453 3196 10 534 3648 10 571
pcb442 2737 11 501 3483 15 414 3989 14 466
gr666 17711 8 789 – – – – – –
mean ratio 1.002535 1.188732
RE A 150 23 8 100 30 7 130 30 6 77
RE A 200 72 7 163 74 8 135 57 6 76
RE A 250 138 9 186 155 7 163 140 6 109
RE A 300 866 6 295 884 6 203 1344 7 211
RE A 350 411 5 252 642 6 147 879 6 150
RE A 400 8456 8 454 6623 7 285 4876 8 296
RE A 450 2107 5 279 1226 4 220 5386 5 215
RE A 500 15330 6 436 13473 6 366 6114 5 237
mean ratio 1.035264 1.291607
mean ratio all 1.016316 1.232048

Table 3: Using no subtours of size 3 vs. using the shortest subtours of size 3 for generation
of subtour constraints before starting the solving process. The parameter p defines the
proportion of used subtour constraints. Mean ratios refer to the arithmetic means over
ratios between the running times of the particular approaches and the running time of
the BasicIntegerTSP (corresponding to p = 0). “sec.” is the time in seconds, “#i.” the
number of iterations and “#c.” the number of subtour elimination constraints added to
the ILP before starting the last iteration. The entries “–” by TSPLIB instances cannot
be computed with 16 GB RAM.

are able to generate during one iteration, but to make a proper selection. We again used
our computational tests in order to identify two general properties which seem to point
to such “suitable” subtour inequalities.

• Sort all obtained subtours with respect to their cardinality, chose the smallest ones
and added the corresponding subtour inequalities into the model.

• Sort all obtained subtours with respect to their length and proceed as above.

8
The corresponding results are summarized in Tables 4 and 5 and it is obvious that
this idea does not speed up our approach as intended. Thus we dropped it from our
considerations.
instance p=1 p = 23 p = 13
sec. #i. #c. sec. #i. #c. sec. #i. #c.
kroA150 19 7 136 34 8 109 69 19 115
kroB150 179 8 148 51 8 135 477 15 134
u159 6 4 49 30 4 52 19 11 56
brg180 44 4 103 27 6 77 59 19 80
kroA200 677 8 237 714 7 171 2846 14 131
kroB200 31 5 121 39 6 98 89 13 77
tsp225 178 9 261 100 14 183 173 34 166
a280 157 11 143 141 12 154 239 27 127
lin318 6885 8 357 7069 12 367 9444 32 392
gr431 2239 9 453 3210 20 522 4924 38 413
pcb442 2737 11 501 1867 18 384 5129 85 386
gr666 17711 8 789 7643 7 505 71594 25 597
mean ratio 1.252892 2.488345
RE A 150 23 8 100 28 9 109 52 17 94
RE A 200 72 7 163 69 8 134 112 23 98
RE A 250 138 9 186 131 10 149 208 20 119
RE A 300 866 6 295 792 10 259 1720 29 293
RE A 350 411 5 252 715 7 232 849 19 177
RE A 400 8456 8 454 129380 8 311 107987 26 299
RE A 450 2107 5 279 1544 7 236 7987 11 238
RE A 500 15330 6 436 18594 8 324 13738 16 308
mean ratio 2.878162 3.354102
mean ratio all 1.903000 2.834648

Table 4: Using all subtours vs. using only the smallest subtours with respect to their
cardinality for generation of subtour constraints. The parameter p defines the propor-
tion of used subtour constraints. Mean ratios refer to the arithmetic means over ratios
between the running times of the particular approaches and the running time of the
BasicIntegerTSP (corresponding to p = 1). “sec.” is the time in seconds, “#i.” the
number of iterations and “#c.” the number of subtour elimination constraints added to
the ILP before starting the last iteration.

3.4 Clustering into subproblems

It can be observed that many subtours have a local context, meaning that a small
subset of vertices separated from the remaining vertices by a reasonably large distance

9
instance p=1 p = 32 p = 31
sec. #i. #c. sec. #i. #c. sec. #i. #c.
kroA150 19 7 136 41 10 131 46 16 90
kroB150 179 8 148 495 7 152 250 16 112
u159 6 4 49 14 5 60 23 12 55
brg180 44 4 103 24 13 86 161 8 78
kroA200 677 8 237 862 6 124 1829 13 132
kroB200 31 5 121 59 7 121 79 11 89
tsp225 178 9 261 112 13 197 197 32 159
a280 157 11 143 94 9 101 212 21 96
lin318 6885 8 357 7688 13 355 9593 36 390
gr431 2239 9 453 6091 15 565 9434 45 530
pcb442 2737 11 501 2365 18 487 5913 70 399
gr666 17711 8 789 14713 10 735 – – –
mean ratio 1.478194 2.434945
RE A 150 23 8 100 24 9 115 45 22 81
RE A 200 72 7 163 60 10 123 113 25 108
RE A 250 138 9 186 138 7 117 209 22 103
RE A 300 866 6 295 1099 10 321 953 23 201
RE A 350 411 5 252 876 7 231 934 16 167
RE A 400 8456 8 454 29625 9 311 301125 27 378
RE A 450 2107 5 279 2926 7 259 4789 14 237
RE A 500 15330 6 436 15786 7 329 37460 16 330
mean ratio 1.524891 6.092589
mean ratio all 1.496873 3.975006

Table 5: Using all subtours vs. using only the smallest subtours with respect to their
length for generation of subtour constraints. The parameter p defines the proportion of
used subtour constraints. Mean ratios refer to the arithmetic means over ratios between
the running times of the particular approaches and the running time of the BasicIn-
tegerTSP (corresponding to p = 1). “sec.” is the time in seconds, “#i.” the number
of iterations and “#c.” the number of subtour elimination constraints added to the
ILP before starting the last iteration. The entries “–” by TSPLIB instances cannot be
computed with 16 GB RAM.

will always be connected by one or more subtours, independently from the size of the
remaining graph (see also Figures 5 and 8 to 19 in the Appendix). Thus, we aim to
identify clusters of vertices and run the BasicIntegerTSP on the induced subgraphs with
the aim of generating within a very small running time the same subtours occurring in
the execution of the approach on the full graph. Furthermore, we can use the optimal
tour from every cluster to generate a corresponding subtour elimination constraint for
the original instance and thus enforce a connection to the remainder of the graph.

10
For our purposes the clustering algorithm should fulfill the following properties:

• clustering quality: The obtained clusters should correspond well to the distance
structure of the given graph, as in a classical geographic clustering.

• running time: Should be low relative to the running time required for the main
part of the algorithm.

• cluster size: If clusters are too large, solving the TSP takes too much time. If
clusters are too small, only few subtour elimination constraints are generated.

Clearly, there is a huge body of literature on clustering algorithms (see e.g. [10])
and selecting one for a given application will never satisfy all our objectives. Our main
restriction was the requirement of using a clustering algorithm which works also if the
vertices are not embeddable in Euclidean space, i.e. only arbitrary edge distances are
given. Simplicity being another goal, we settled for the following approach described in
Algorithm 2:

Input: Complete graph G = (V, E), where |V | = n and |E| = m = n(n−1) 2 , distance
+
function d : E → R0 and parameter c, where 1 ≤ c ≤ n.
Output: Clustering C = {V1 , . . . , Vc }, where V1 ∪ . . . ∪ Vc = V .
1: sort the edges such that de1 ≤ . . . ≤ dem ;
2: define G′ = (V ′ , E ′ ) such that V ′ = V and E ′ = ∅;
3: let i ..= 1;
4: define C ..= {v1 }, . . . , {vn } ;

5: while |C| > c do

6: set E ′ ..= E ′ ∪ {ei };
7: set C ..= {V1 , . . . , V|C| }, where V1 , . . . , V|C| are the connected components of graph
G′ ;
8: end while

Algorithm 2: Clustering algorithm.

First, we fix the number of clusters c with 1 ≤ c ≤ n and sort the edges in increasing
order of distances (see line 1). Then we start with the empty graph G′ = (V ′ , E ′ ) (line 2)
containing only isolated vertices (i.e. n clusters) and add iteratively edges in increasing
order of distances until the desired number of clusters c is reached (see lines 5 and 6).
In each iteration the current clustering is implied by the connected components of the
current graph (see line 7). We denote this clustering approach by C | c. Note that this
clustering algorithm does not make any assumptions about the underlying TSP instance
and does not exploit any structural properties of the Metric TSP or the Euclidean TSP.
It was observed in our computational experiments that the performance of the TSP
algorithm is not very sensitive to small changes of the cluster number c and thus a rough
estimation of c is sufficient. The behavior of the running time as a function of c can be
found for particular test instances in Figure 21, see Section 4.2 for further discussion.

11
3.5 Restricted clustering
Although the clustering algorithm (see Algorithm 2) decreases the computational time
of the whole solution process for some test instances, we observed a certain shortcoming.
There may easily occur clusters consisting of isolated points or containing only two ver-
tices. Clearly, these clusters do not contribute any subtour on their own. Moreover, the
degree constraints (2) guarantee that each such vertex is connected to the remainder of
the graph in any case. The connection of these vertices to some “neighboring” cluster
enforced in BasicIntegerTSP implies that the clustering yields different subtours for
these neighbors and not the violated subtour elimination constraints arising in BasicIn-
tegerTSP.
To avoid this situation, we want to impose a minimum cluster size of 3. An easy way
to do so is as follows: After reaching the c clusters, continue to add edges in increasing
order of distances (as before), but add an edge only, if it is incident to one of the vertices
in a connected component (i.e. cluster) of size one or two. This means basically that we
simply merge these small clusters to their nearest neighbor with respect to the actual
clustering. Note that this is a step-by-step process and it can happen that two clusters
of size 1 merge first before merging the resulting pair to its nearest neighboring cluster.
The resulting restricted clustering approach will be denoted by RC3 | c.
Against our expectations, the computational experiments (see Section 4) show that
this approach often impacts the algorithm in the opposite way (see also Figure 21 and
Table 9 in the Appendix) if compared for the same original cluster size c.
Surprisingly, we could observe an interesting behavior if c ≈ n. In this case, the
main clustering algorithm (see Algorithm 2) has almost no effect, but the “post-phase”
which enforces the minimum cluster size yields a different clustering on its own. This
variant often beats the previous standard clustering algorithm with c ≪ n (see Table 9
in the Appendix). Note that we cannot fix the actual number of clusters c′ in this
case. But our computational results show that c′ ≈ n5 usually holds if the points are
distributed relatively uniformly in the Euclidean plane and if the distances correspond
to their relative Euclidean distances (see Figure 20 in the Appendix).

3.6 Hierarchical clustering

It was pointed out in Section 3.4 that the number of clusters c is chosen as an input
parameter. The computational experiments in Subsection 4.2 give some indication on
the behavior of Algorithm 2 for different values of c, but fail to provide a clear guideline
for the selection of c. Moreover, from graphical inspection of test instances, we got the
impression that a larger number of relevant subtour elimination constraints might be
obtained by considering more clusters of moderate size. In the following we present an
idea that takes both of these aspects into account.
In our hierarchical clustering process denoted by HC we do not set a cluster number c,
but let the clustering algorithm continue until all vertices are connected (this corresponds
to c = 1). The resulting clustering process can be represented by a binary clustering
tree which is constructed in a bottom-up way. The leaves of the tree represent isolated

12
vertices, i.e. the n trivial clusters given at the beginning of the clustering algorithm.
Whenever two clusters are merged by the addition of an edge, the two corresponding
tree vertices are connected to a new common parent vertex in the tree representing
the new cluster. At the end of this process we reach the root of the clustering tree
corresponding to the complete vertex set. An example of such a clustering tree is shown
in Figures 1 and 2.

v3

v1 v2 v4 v5

Figure 1: Example illustrating the hierarchical clustering: Vertices of the TSP instance.
Distances between every two vertices correspond to their respective Euclidean distances
in this example.

{v1 , v2 , v3 , v4 , v5 }

{v1 , v2 , v3 }

{v1 , v2 } {v4 , v5 }

{v1 } {v2 } {v3 } {v4 } {v5 }

Figure 2: Example illustrating the hierarchical clustering: Clustering tree.

Now, we go through the tree in a bottom-up fashion from the leaves to the root.
In each tree vertex we solve the TSP for the associated cluster, after both of its child
vertices were resolved. The crucial aspect of our procedure is the following: All subtour
elimination constraints generated during such a TSP solution for a certain cluster are
propagated and added to the ILP model used for solving the TSP instance of its parent
cluster. Obviously, at the root vertex the full TSP is solved.
The advantage of this strategy is the step-by-step construction of the violated subtour

13
elimination constraints. A disadvantage is that many constraints can make sense in the
local context but not in the global one and thus too many constraints could be generated
in this way. Naturally, one pays for the additional subtour elimination constraints by an
increase in computation time required to solve a large number of – mostly small – TSP
instances. To avoid the solution of TSPs of the same order of magnitude as the original
instance, it makes sense to impose an upper bound u on the maximum cluster size. This
means that the clustering tree is partitioned into several subtrees by removing all tree
vertices corresponding to clusters of size greater than u. After resolving all these subtrees
we collect all generated subtour elimination constraints and add them to the ILP model
for the originally given TSP. This approach will be denoted as HC | u. Computational
experiments with various choices of u indicated that u = 4 logn n would be a good upper
2
bound.
Let us take a closer look at the problem of including too many subtour elimination
constraints which are redundant in the global graph context. Of course the theoretical
“best” way would be to check which of the propagated subtour elimination constraints
were not used during the runs of the ILP solver and drop them. To do this, it would be
necessary to get this information from the ILP solver which often is not possible.
However, we can try to approximately identify subtours which are not only locally
relevant in the following way: All subtour elimination constraints generated in a certain
tree vertex, i.e. for a certain cluster, are marked as considered subtour elimination con-
straints. Then we solve the TSP for the cluster of its parent vertex in the tree without
using the subtours marked as considered. If we generate such a considered subtour again
during the solution of the parent vertex, we take this as an indicator of global significance
and add the constraint permanently for all following supersets of this cluster. If we set
the upper bound u, we take also all subtour elimination constraints found in the biggest
solved clusters. This approach will be denoted as HCD | u.
Of course, it is only a heuristic rule and one can easily find examples, where this
prediction on a subtour’s relevance fails, but our experiments indicate that HCD |
4n/ log2 n is the best approach we considered. A comparison with other hierarchical
clustering methods for all test instances can be found in Table 8 in the Appendix. It
can be seen that without an upper bound we are often not able to find the solution at
all (under time and memory constraints we made on the computational experiments).
In the third and fourth column we can see a comparison between approaches both using
the upper bound u = 4 logn n where the former collects all detected subtour elimination
2
constraints and the latter allows to drop those which seem to be relevant only in a local
context. Both these methods beat BasicIntegerTSP (for the comparison of this approach
with other presented algorithms see the computational experiments in Section 4).

4 Computational experiments
In the following the computational experiments and their results will be discussed.

14
4.1 Setup of the computational experiments
All tests were run on an Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz with 16 GB RAM
under Linux1 and all programs were implemented in C++2 by using the SCIP MIP-
solver [1] together with CPLEX as LP-solver3 . It has often been discussed in the lit-
erature (see e.g. [13]) and in personal communications that ILP-solvers are relatively
unrobust and often show high variations in their running time performance, even if the
same instance is repeatedly run on the same hardware and same software environment.
Our first test runs also exhibited deviations up to a factor of 2 when identical tests
were repeated. Thus we took special care to guarantee the relative reproducibility of the
computational experiments: No additional swap memory was made available during the
tests, only one thread was used and no other parallel user processes were allowed. This
leads to a high degree of reproducibility in our experiments. However, this issue makes
a comparison to other simple approaches, which were tested on other computers under
other hardware and software conditions, extremely difficult.
We used two groups of test instances: The first group is taken from the well-known
TSPLIB95 [17], which contains the established benchmarks for TSP and related prob-
lems. From the collection of instances we chose all those with (i) at least 150 and at most
1000 vertices and (ii) which could be solved in at most 12 hours by our BasicIntegerTSP.
It turned out that 25 instances of the TSPLIB95 fall into this category (see Table 9),
the largest having 783 vertices.
We also observed some drawbacks of these instances: Most of them (23 of 25) are
defined as point sets in the Euclidean plane with distances corresponding to the Euclidean
metric or as a set of geographical cities, i.e. points on a sphere. Moreover, they often
contain substructures like meshes or sets of colinear points and finally, since all distances
are rounded to the nearest integer, there are many instances which have multiple optimal
solutions. These instances are relatively unstable with respect to solution time, number
of iterations, and – important for our approach – cardinality of the set of violated subtour
elimination constraints. For our approach instances with a mesh geometry (e.g. ts225
from TSPLIB95) were especially prone to unstable behavior, such as widely varying
running times for minor changes in the parameter setting. This seems to be due to the
fact that these instances contain many 2-matchings with the same objective function
value as illustrated in the following example: Consider a 3 × (2n + 2) mesh graph (see
Figure 3, left graph). It has 2n optimal TSP tours (see [21]). If we fix a subtour
on the first 6 vertices,
we obviously have 2n−1 optimal TSP tours on the remaining
3 × (2n + 2) − 2 vertices (see Figure 3, right graph) and together with the fixed subtour
we have 2n−1 2-matchings having the same objective value as an optimal TSP tour on
the original graph. Thus the search process for a feasible TSP tour can vary widely.
In order to provide further comparisons, we also defined a set of instances based
1
precise version: Linux 3.8.0-29-generic #42˜precise1-Ubuntu SMP x86 64 x86 64 x86 64 GNU/Linux
2
precise compiler version: gcc version 4.6.3
3
precise version: SCIP version 3.0.1 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver:
CPLEX 12.4.0.0] [GitHash: 9ee94b7] Copyright (c) 2002-2013 Konrad-Zuse-Zentrum für Informations-
technik Berlin (ZIB)

15
3 3

2n + 2 (2n + 2) − 2

Figure 3: Example illustrating the behavior of our approaches by instances based on

graphs containing mesh substructures. Distances between every two vertices correspond
to their respective Euclidean distances in this example.

on random Euclidean graphs: In a unit square [0, 1]2 we chose n uniformly distributed
points and defined the distance between every two vertices as their respective Euclidean
distance4 . These random Euclidean instances eliminate the potential influence of sub-
structures and always have only one unique optimal solution in all stages of the solving
process. We created 40 such instances named RE X n where n ∈ {150, 200, 250, . . . , 500}
indicates the number of vertices and X ∈ {A, B, C, D, E}.
The running times of our test instances, most of them containing between 150 and
500 vertices, were often within several hours. Since we tested many different variants
and configurations of our approach, we selected a subset of these test instances to get
faster answers for determining the best algorithm settings for use in the final tests. This
subset contains 12 (of the 25) TSPLIB instances and one random instances for every
number of vertices n (see e.g. Table 1.)
All our running time tables report the name of the instance, the running time (sec.)
in wall-clock seconds (rounded down to nearest integers), the number of iterations (#i.),
i.e. the number of calls to the ILP-solver in the main part of our algorithm (without the
TSP solutions for the clusters) and the number of subtour elimination constraints (#c.)
added to the ILP model in the last iteration, i.e. the number of constraints of the model
which yielded an optimal TSP solution. We often compare two columns of a table by
taking the mean ratio, i.e. computing the quotient between the running times on the
same instance and taking the arithmetic mean of these quotients.

4.2 Computational details for selected examples

Let us now take a closer look at two instances in detail. While this serves only as an
illustration, we studied lots of these special case scenarios visually during the development
of the clustering approach to gain a better insight into the structure of subtours generated
by BasicIntegerTSP.
We selected instances kroB150 and u159 whose vertices are depicted in Figures 6
and 7 in the Appendix. Both instances consist of points in the Euclidean plane and the
4
We represented all distances as integers by scaling with 214 and rounding to the nearest integer.

16
distances between every two vertices correspond to their respective Euclidean distances,
however, they represent two very different instance types: The instance kroB150 consists
of relatively uniformly distributed points, the instance u159 is more structured and it
contains e.g. mesh substructures which are the worst setting for our algorithm (recall
Subsection 4.1).
Figure 21 in the Appendix illustrates the behavior of the running time t in seconds as a
function of the parameter c for the instances kroB150 and u159. The full lines correspond
to standard clustering approach C | c described in Section 3.4 (see Algorithm 2), while
the dashed line corresponds to the restricted clustering RC3 | c of Section 3.5 with
minimum cluster size 3. The standard BasicIntegerTSP without clustering arises for
c = 1.
Instance kroB150 consists of relatively uniformly distributed points in the Euclidean
plane, but has a specific property: By using Algorithm 2 we can observe the occurrence
of two main components also for relatively small coefficient c (already for c = 6). This
behavior is rather atypical for random Euclidean graphs, cf. [16, ch. 13], but it provides
an advantage for our approach since we do not have to solve cluster instances of the same
order of magnitude as the original graph but have several clusters of moderate size also
for small cluster numbers c.
Considering the standard clustering approach (Algorithm 2) in Figure 21, upper
graph, it can be seen that only a small improvement occurs for c between 2 and 5.
Looking at the corresponding clusterings in detail, it turns out in these cases that there
exists only one “giant connected component” and all other clusters have size 1. This
structure also implies that for the restricted clustering these isolated vertices are merged
with the giant component and the effect of clustering is lost completely. For larger cluster
numbers c, a considerable speedup is obtained, with some variation, but more or less in
the same range for almost all values of c ≥ 6 (in fact, the giant component splits in these
cases). Moreover, the restricted clustering performs roughly as good as the standard
clustering for c ≥ 6.
Instance u159 is much more structured and has many colinear vertices. Here, we
can observe a different behavior. While the standard clustering is actually beaten by
BasicIntegerTSP for smaller cluster numbers and has a more or less similar performance
for larger cluster numbers, the restricted clustering is almost consistently better than
the other two approaches. For c between 2 and 10 there exists a large component
containing many mesh substructures which consumes as much computation time as the
whole instance.
These two instances give some indication of how to characterize “good” instances for
our algorithm: They should
• consist of more clearly separated clusters and
• not contain mesh substructures and colinear vertices.

4.3 General computational results

A summary of the computational results for BasicIntegerTSP and the most promising
variants of clustering based subtour generations can be found in Table 9. For random

17
Euclidean instances we report only the mean values of all five instances of the same
size. It turns out that HCD | 4 logn n , i.e. the hierarchical clustering approach combined
2
with dropping subtour elimination constraints and fixing them only if they are generated
again in the subsequent iteration and with the upper bound on the maximum cluster
size u = 4 logn n , gives the best overall performance. A different behavior can be observed
2
for instances taken from the TSPLIB and for random Euclidean instances. On the
TSPLIB instances this algorithm HCD | 4 logn n is on average about 20% faster than pure
2
BasicIntegerTSP and beats the other clustering based approaches for most instances. In
those cases, where it is not the best choice, it is usually not far behind.
As already mentioned, best results are obtained with HCD | 4 logn n for instances with
2
a strong cluster structure and without mesh substructures (e.g. pr299). For instances
with mesh substructures it is difficult to find an optimal 2-matching which is also a TSP
tour. For random Euclidean instances the results are less clear but approaches with fixed
number of clusters seem to be better then the hierarchical ones.
It was a main goal of this study to find a large number of “good” subtour elimin-
ation constraints, i.e. subtours that are present in the last iteration of the ILP-model
of BasicIntegerTSP. Therefore, we show the potentials and limitations of our approach
in reaching this goal. In particular, we will report the relation between the set S1 con-
sisting of all subtours generated by running a hierarchical clustering algorithm with an
upper bound u (set as in the computational tests to u = 4 logn n ) before solving the
2
original problem (i.e. the root vertex) and the set S2 containing only the subtour elim-
ination constraints included in the final ILP model of BasicIntegerTSP. We tested the
hierarchical clustering with and without the dropping of non-repeated subtours.
There are two aspects we want to describe: At first, we want to check whether
S1 contains a relevant proportion of “useful” subtour contraints, i.e. constraints also
included in S2 , or whether S1 contains “mostly useless” subtours. Therefore, we report
the proportion of used subtours defined as
|S1 ∩ S2 |
pused ..= . (7)
|S1 |
Secondly, we want to find out to what extend it is possible to find the “right” subtours
by our approach. Hence, we define the proportion of covered subtours defined as
|S1 ∩ S2 |
pcov ..= . (8)
|S2 |
The values of pused and pcov are given in Table 6. It can be seen that empirically there
is the chance to find about 26–31% (pcov ) of all required violated subtour elimination
constraints. If subtour elimination constraints are allowed to be dropped, we are able to
find fewer such constraints, but our choice has a better quality (pcov is smaller, but pused
is larger), i.e. the solver does not have to work with a large number of constraints which
only slow down the solving process and are not necessary to reach an optimal solution.
Furthermore, we can observe a relative big difference between the values of the pro-
portion of used subtour elimination constraints (pused ) for the TSPLIB instances and for
random Euclidean instances if the dropping of redundant constraints is allowed.

18
instance HC | 4 logn n HCD | 4 logn n
2 2
pused pcov pused pcov
kroA150 0.262712 0.455882 0.476190 0.367647
kroB150 0.222222 0.351351 0.396040 0.270270
u159 0.085271 0.448980 0.153226 0.387755
brg180 0.133929 0.145631 0.714286 0.145631
kroA200 0.209713 0.324895 0.450704 0.270042
kroB200 0.206612 0.413223 0.423423 0.388430
tsp225 0.134752 0.218391 0.297143 0.199234
a280 0.064935 0.314685 0.161943 0.279720
lin318 0.234589 0.383754 0.440273 0.361345
gr431 0.073701 0.209713 0.221053 0.185430
pcb442 0.056759 0.151697 0.133117 0.163673
gr666 0.076048 0.271229 0.220379 0.235741
mean 0.146770 0.307453 0.340648 0.271243
RE A 150 0.179191 0.310000 0.289157 0.240000
RE A 200 0.122642 0.239264 0.212329 0.190184
RE A 250 0.120773 0.268817 0.172727 0.204301
RE A 300 0.191235 0.325424 0.331915 0.264407
RE A 350 0.151274 0.376984 0.285714 0.333333
RE A 400 0.170455 0.297357 0.254157 0.235683
RE A 450 0.148148 0.415771 0.311178 0.369176
RE A 500 0.165485 0.321101 0.276596 0.268349
mean 0.156150 0.319340 0.266722 0.263179
mean of all 0.150522 0.312207 0.311078 0.268018

Table 6: Proportion of used and proportion of covered subtours for our hierarchical
clustering approaches with the upper bound u = 4 logn n which (i) does not allow (HC |
2
4 logn n ) and which (ii) does allow (HCD | 4 logn n ) to drop the unused subtour elimination
2 2
constraints.

4.4 Adding a starting heuristic

Of course, there are many possibilities of adding improvements to our basic approach.
Lower bounds and heuristics can be introduced, branching rules can be specified, or
cutting planes can be generated. We did not pursue these possibilities since we want to
focus on the simplicity of the approach. Moreover, we wanted to take the ILP solver as
a “black box” and not interfere with its execution.
Just as an example which immediately comes to mind, we added a starting heuristic
to give a reasonably good TSP solution as a starting solution to the ILP solver. We
used the improved version of the classical Lin-Kernighan heuristic in the code written
by Helsgaun [9]. The computational results reported in Table 7 show that a considerable
speedup (roughly a factor of 3, but also much more) can be obtained in this way.

19
instance without starting with starting
heuristic heuristic
sec. #i. #c. sec. #i. #c.
kroA150 19 7 136 16 10 34
kroB150 179 8 148 17 8 104
u159 6 4 49 4 5 40
brg180 44 4 103 0 2 15
kroA200 677 8 237 42 8 135
kroB200 31 5 121 28 6 124
tsp225 178 9 261 73 13 176
a280 157 11 143 32 8 58
lin318 6885 8 357 4941 8 259
gr431 2239 9 453 838 10 318
pcb442 2737 11 501 447 18 207
gr666 17711 8 789 13225 11 485
mean ratio (sec.) 0.432074
RE A 150 23 8 100 14 11 65
RE A 200 72 7 163 38 11 99
RE A 250 138 9 186 63 9 124
RE A 300 866 6 295 146 8 173
RE A 350 411 5 252 126 6 151
RE A 400 8456 8 454 1274 6 251
RE A 450 2107 5 279 482 7 197
RE A 500 15330 6 436 1997 9 241
mean ratio (sec.) 0.322231
mean ratio all 0.388137

Table 7: Results for BasicIntegerTSP used without / with the Lin-Kernighan heuristic
for generating an initial solution. Mean ratios refer to the arithmetic means over ratios
between the running times of the approaches. “sec.” is the time in seconds, “#i.” the
number of iterations and “#c.” the number of subtour elimination constraints added to
the ILP before starting the last iteration.

5 Some theoretical results and further empirical observa-

tions
Although our work mainly aims at computational experiments, we also tried to analyze
BasicIntegerTSP from a theoretical point of view. In particular we studied the expected
behavior on random Euclidean instances and tried to characterize the expected cardin-
ality of the minimal set of required subtours S ∗ as defined in Section 2. It is well known
that no polynomially bounded representation of the TSP polytope can be found and there
also exist instances based on a mesh-structure for which E [|S ∗ |] has exponential size, but

20
the question for the expected size of |S ∗ | for random Euclidean instances and thus for
the expected number of iterations of our solution algorithm remains an interesting open
problem.
We started with extensive computational tests, some of them presented in Figures 22
and 23 in the Appendix, to gain empirical evidence on this aspect. The upper graph
in Figure 22 illustrates the mean number of iterations needed by BasicIntegerTSP to
reach optimality for different numbers of vertices n (we evaluated 100 random Euclidean
instances for every value n). The lower graph of Figure 22 shows the mean length of the
optimal TSP tour and of the optimal 2-matching (i.e. the objective value after solving
the ILP in the first iteration) by using the same setting.
It was proven back in 1959 that the expected length of an optimal TSP tour is
√
asymptotically β n, where β is a constant [3]. This approach was later generalized for
other settings and other properties of the square root asymptotic were identified [19, 22].
We used these properties to prove the square root asymptotic also for the 2-matching
problem (cf. Figure 22, lower graph, dashed).
We need some definitions definitions, lemmas and theorems originally introduced
by [19] and summarized by [22] first in order to prove this result.
Definition 1 ((2-)matching functional and boundary (2-)matching functional).
Let F ..= F(dim) denote the finite subsets of Rdim and let let R ..= R(dim) denote
the dim-dimensional rectangles of Rdim .
Furthermore, let F ∈ F be a point set in Rdim and let R ∈ R be a dim-dimensional
rectangle in Rdim where dim ≥ 2.
And finally, let d : R × R → R+

0 be a metric and let G = G(F, R) = V (G), E(G) be
a complete graph with the vertex set V (G) = F ∩ R and with the distances d(e) between
every two vertices u, v ∈ V (G) where e = (u, v).
Then we will denote
M (F, R) ..= M (F ∩ R) ..= min {OV (G, m)|m is a matching in G} (9)
m
the matching functional.
Furthermore, we will denote
 
 ( )
 X 
MB (F, R) ..= MB (F ∩ R) ..= min M (F, R), inf M (Fi ∪ {ai , bi }, R)
 (Fi )i≥1 ∈F 
 i 
{ai ,bi }i≥1 ∈∂R

(10)
the boundary matching functional. F stays for the set of all partitions (Fi )i≥1 of F
and ∂R stays for the set of all sequences of pairs of points {ai , bi }i≥1 belonging to the
boundary of the rectangle R denoted by ∂R. Additionally, we set d(a, b) = 0 for all
a, b ∈ ∂R.
Similarly, we define the 2-matching functional L and the boundary 2-matching func-
tional LB .
L(F, R) ..= L(F ∩ R) ..= min {OV (G, x)|x is a 2-matching in G)} (11)
x

21
 
 ( )
 X 
LB (F, R) ..= LB (F ∩R) ..= min L(F, R), inf L(Fi ∪ {ai , bi }, R) (12)
 (Fi )i≥1 ∈F 
 i 
{ai ,bi }i≥1 ∈∂R

If it is obvious which rectangle R is considered, we also write M (F ) for M (F, R),

MB (F ) for MB (F, R), L(F ) for L(F, R) and LB (F ) for LB (F, R).
Finally, we define L(F, R) = LB (F, R) = 0 if |F ∩ R| < 3
Definition 2 (simple subadditivity, geometric subadditivity and superadditivity). Let
R be a rectangle defined as [0, t]dim for some positive constant t ∈ R+ partitioned into
two rectangles R1 and R1 (R1 ∪ R2 = R). Furthermore, let F, G be finite sets in [0, t]dim
and let P (F, R) : F × R → R+0 be a function.
The function P is simple subadditive if the following inequality is satisfied:

P (F ∪ G, R) ≤ P (F, R) + P (G, R) + C1 t (13)

where C1 ..= C1 (dim) is a finite constant.

If

P (F, R) ≤ P (F, R1 ) + P (F, R2 ) + C2 diam(R) (14)

is fulfilled, we will call the function P geometric subadditive. diam(R) denotes the dia-
meter of the rectangle R and C2 ..= C2 (dim) is a finite constant.
If

P (F, R) ≥ P (F, R1 ) + P (F, R2 ) (15)

is satisfied, we will call the function P superadditive.

Definition 3 (subadditive and superadditive Euclidean functional). Let
P (F, R) : F × R → R+
0 be a function satisfying

∀R ∈ R, P (∅, R) = 0 , (16)

∀y ∈ Rdim , R ∈ R, F ⊂ R : P (F, R) = P (F + y, R + y) , (17)

∀α > 0, R ∈ R, F ⊂ R : P (αF, αR) = αP (F, R) . (18)

Then we will say that P is translation invariant (condition (17)) and homogeneous (con-
dition (18)).
We will call P an Euclidean functional.
Let R be a rectangle defined as [0, t]dim for some positive constant t ∈ R+ partitioned
into two rectangles R1 and R1 (R1 ∪ R2 = R). If P satisfies the geometric subadditivity
(14), we will say that P is a subadditive Euclidean functional. If P is supperadditive
(15), we will say that P is a superadditive Euclidean functional.

22
Lemma 4 ([22], originally [19]). The matching functional M and the boundary matching
functional MB are subadditive Euclidean functionals.

Proof. See [22].

Lemma 5 (growth bounds – [22], originally [19]). Let P be a subadditive Euclidean

functional. Then there exists a finite constant C4 ..= C4 (dim) such that for all dim-
dimensional rectangles of Rdim and for all F ⊂ R we have
dim−1
P (F, R) ≤ C4 diam(R)|F | dim . (19)

Proof. See [22].

Definition 6 (smoothness). An Euclidean functional P is smooth if there is a finite

constant C3 ..= C3 (dim) such that for all sets F, G ∈ [0, 1]dim we have

P (F ∪ G) − P (F ) ≤ C3 |G| d−1

d . (20)

Definition 7 (pointwise closeness). Say that Euclidean functionals P and PB are point-
wise close if for all subsets F ⊂ [0, 1]dim we have

P F, [0, 1]dim − PB F, [0, 1]dim = o |F | d
d−1
. (21)

Definition 8 (complete convergence). A sequence of random variables Xn , n ≥ 1,

converges completely (c.c.) to a constant C if and only if for all ε > 0 we have
∞
X
P [|Xn − C| > ε] < ∞ . (22)
n=1

Theorem 9 (basic limit theorem for Euclidean functionals – [22], originally [19]). Let
Xi , 1 ≤ i ≤ n, be independent and identically distributed random variables with values
in the unit dim-dimensional rectangle [0, 1]dim , dim ≥ 2.
If PB is a smooth superadditive Euclidean functional on Rdim , dim ≥ 2, then
PB (X1 , X2 , . . . , Xn )
lim = α(PB , dim) c.c. , (23)
n→∞ dim−1
n dim
where α(PB , dim) is a positive constant.
If P is an Euclidean functional on Rdim , dim ≥ 2, which is pointwise close to PB ,
then
P (X1 , X2 , . . . , Xn )
lim = α(PB , dim) c.c. (24)
n→∞ dim−1
n dim
Proof. See [22].

23
We can prove our result now.

Lemma 10. The 2-matching functional L and the boundary 2-matching functional LB
fulfill the conditions of Theorem 9.

Proof. The proof is a modification of similar proofs for other combinatorial optimization
problems contained in [22].

(1) First, we show that the boundary 2-matching functional LB is a superadditive

Euclidean functional. Equalities (16), (17) and (18) are fulfilled obviously.
Let us now show the superadditivity. We can distinguish 2 cases in general:

(a) Either the solution over the whole rectangle R does not cross the boundary
between the rectangles R1 and R2 at all or
(b) at least one subtour crosses the boundary between the rectangles R1 and R2 .

LB (F, R) = LB (F, R1 ) + LB (F, R2 ) obviously holds in the first case.

Let us now consider a subtour crossing the boundary between the rectangles R1
and R2 (for an example see Figure 4). W.l.o.g. we can assume that the boundary
is crossed between the points v1 and v2 , and v3 and v4 and that the crossing points
are x and y respectively. Furthermore, w.l.o.g. we can assume that v1 , v3 ∈ R1 .
Then the new subtour, containing the vertices v1 and v3 , and lying in the rectangle
R1 , consists of the following parts:

• the same path between the vertices v1 and v3 belonging to the rectangle R1
as in the whole rectangle R,
• the orthogonal connections between the vertices v1 and v3 and the boundary
between the both rectangles, and finally
• a piece of this boundary (see also Figure 4).

We have to choose the vertices a and b on this boundary in such a way that
a = arg min {d(v1 , α)} and b = arg min {d(v3 , β)} hold in order to achieve the
α∈∂R1 ∩∂R2 β∈∂R1 ∩∂R2
minimality. Due to this choice of the vertices a and b we can write d(v1 , a) ≤ d(v1 , x)
and d(v3 , b) ≤ d(v3 , y) and since d(a, b) = 0 and the remaining part of the subtour
belonging to the rectangle R1 yield the same contribution to the objective value,
we can claim that the contribution of this new subtour to the objective value is
smaller or equal to the contribution of the part of the original subtour lying in the
rectangle R1 . The same argument can be used for the second rectangle R2 and for
all other subtours crossing the boundary between the rectangles R1 and R2 .

(2) Next, we check some properties of the 2-matching functional L. Equalities (16),
(17) and (18) obviously hold.
Further, it is easy to see that the 2-matching functional L fulfill the geometric
subadditivity. Since we minimize, the minimum weighted 2-matching over the

24
v1
a
x
v2 v10
a′ v11

R1 R2 v9
v8
v7
b
v3
v6 y

v5 b′ v4

Figure 4: Example illustrating the superadditivity of the boundary 2-matching functional

LB .

whole rectangle R can have only a smaller objective value than the sum of the
objective values for the rectangles R1 and R2 taken separately.
As a next step, we have to prove the pointwise closeness of the 2-matching func-
tional L to the boundary 2-matching functional LB .

First, note that LB F, [0, 1]dim ≤ L F, [0, 1]dim always hold (see (12)). Thus
it suffices to show

L F, [0, 1]dim ≤ LB F, [0, 1]dim + C7 |F | d

d−1
(25)

where C7 ..= C7 (dim) is a finite constant.

Let F ∗ ⊆ F be the set of vertices which are connected with the boundary ∂[0, 1]dim
by a path. Now, we remove all edges incident with the vertices contained in the
vertex set F ∗ . If |F ∗ | < 3, we can just put these vertices to an arbitrary subtour (if
such a subtour exists) and get the √ above inequality (the increase of the objective
value can be easily bounded by 4 dim in this case). If |F ∗ | ≥ 3, we can construct
an minimum weighted 2-matching on this vertex set and obtain

L F, [0, 1]dim ≤ LB F, [0, 1]dim + L F ∗ , [0, 1]dim .

(26)

And since |F ∗ | ≤ |F |, we can use Lemma 5 and get inequality (25).

(3) Finally, we prove the smoothness of the boundary 2-matching functional LB . We

will show the simple and geometric subadditivity of this functional first in order to
be ably to show the smoothness.
Let F and G be finite sets in [0, t]dim . If the minimum weighted 2-matching for
the vertex set F ∪ G equals to the minimum weighted 2-matchings for the vertex

25
sets F and G joined together, inequality (13) holds with equality. Since we can
always construct such a solution for the vertex set F ∪ G, the objective value can
be only smaller in the other case (note that we minimize it).
Let us now prove the geometric subadditivity in order to fulfill the conditions of
Lemma 5. We know that the 2-matching functional L is geometric subadditive.
Now, it is easy to see that the proof of inequality (25) can be easily modified in
order to obtain
d−1
L(F, R) ≤ LB (F, R) + C7 diam(R)|F | d (27)

for an arbitrary dim-dimensional rectangle. Since LB (F, R) ≤ L(F, R), we obtain

inequality (14) immediately.
Using the simple subadditivity and Lemma 5 we get for all finite sets F, G ⊂
[0, 1]dim
√ dim−1
LB (F ∪ G) ≤ LB (F ) + C1 + C4 dim |G| dim
(28)
dim−1
≤ LB (F ) + C5 |G| dim

where C5 ..= C5 (dim) denotes a finite constant. This completes this part of the
proof if LB (F ∪ G) − LB (F ) ≥ 0. Hence we just need to show the following
inequality
dim−1
LB (F ∪ G) ≥ LB (F ) − C6 |G| dim (29)

for some finite constant C6 ..= C6 (dim).

Consider the global minimum weighted 2-matching on the vertex set G ∪ F and
remove all edges from all subtours incident with a vertex g ∈ G. This yield at
most |G| paths of a length of at least 1 containing only vertices from the vertex
set F and some isolated points F ′ ⊆ F . Let F ∗ denote the set of all endpoints
of those paths. Clearly, we have |F ′ | ≤ |G| and |F ∗ | ≤ 2|G|. Consider now the
boundary matching functional MB (F ∗ ) and the corresponding matching m. This
matching together with the disconnected paths and with parts of the boundary
of the dim-dimensional rectangle [0, 1]dim yield a set of subtours {F̃i }N
i=1 for some
particular positive integer N . Furthermore, we can construct an minimum weighted
2-matching on the vertex set F ′ and get a feasible minimum weighted 2-matching.
We can write

LB (F ) ≤ LB (F ∪ G) + MB (F ∗ ) + LB (F ′ ) . (30)

By using Lemmas 4 and 5 we obtain

!
dim−1 dim−1
LB (F ) ≤ LB (F ∪ G) + C6 |F | dim + |F ′ | dim
∗
. (31)

26
And since |F ′ | ≤ |G| ≤ 2|G| and |F ∗ | ≤ 2|G|, we get
!
dim−1 dim−1
LB (F ) ≤ LB (F ∪ G) + C6 (2|G|) dim + 2 (|G|) dim
(32)
dim−1
≤ LB (F ∪ G) + C6 (|G|) dim .

This is exactly inequality (29).

Theorem 11. Let G = (V, E) be a random Euclidean graph with n = |V | vertices and
let d : E → R+0 be the Euclidean distance function. Furthermore, let M2 (G, d) be the
length of an optimal 2-matching. Then

M2 (G, d)
lim √ = α c.c., where α > 0. (33)
n→∞ n

Proof. The theorem immediately follows from Theorem 9 and Lemma 10.

Based on these results the following idea might lead to a proof that the expected
cardinality S ∗ is polynomially bounded: After the first iteration of the algorithm we
have a solution possibly consisting of several separate subtours of total asymptotic length
√ √
α n = α1 n. If there are subtours, we add subtour elimination constraints (in fact at
most ⌊ n3 ⌋), resolve the enlarged ILP and get another solution whose asymptotic length
√
is α2 n. By proving that the expected length of the sequence α = α1 , . . . , α#i = β is
polynomially bounded in n, one would obtain that also E [|S ∗ |] is polynomially bounded
since only polynomially many subtours are added in each iteration. Our intuition and
computational tests illustrated in Figure 22, upper graph, indicate that the length of
√
this sequence could be proportional to n as well. Unfortunately, we could not find the
suitable techniques to show this step.
A different approach is illustrated by Figure 23, where we examine the mean number
of subtours contained in every iteration. In particular, we chose n = 60, generated 100000
random Euclidean instances and sorted them by the number of iterations #i. required by
BasicIntegerTSP. The most frequent number of ILP solver runs was 7 (dotted line), but
we summarize the results for 5 (full line), 6 (dashed), 8 (loosely dashed) and 9 (loosely
dotted) necessary runs in this figure as well. For every iteration of every class (with
respect to the number of involved ILP runs) we compute the mean number of subtours
contained in the respective solutions. As can be expected these numbers of subtours are
decreasing (in average) over the number of iterations. To allow a better comparison of
this behavior for different numbers of iterations we scaled the iteration numbers into the
interval [0, 10] (horizontal axis of Figure 23). It can be seen that the average number
of subtours contained in an optimal 2-matching (first iteration) is about 9.2 while in
the last iteration we trivially have only one tour. Between these endpoints we can first
observe a mostly convex behavior, only in the last step before reaching the optimal TSP

27
tour a sudden drop occurs. It would be interesting to derive an asymptotic description
of these curves. An intuitive guess would point to an exponential function, but so far we
could not find a theoretical justification of this claim.

6 Conclusions
In this paper we provide a “test of concept” of a very simple approach to solve TSP
instances of medium size to optimality by exploiting the power of current ILP solvers.
The approach consists of iteratively solving ILP models with relaxed subtour elimina-
tion constraints to integer optimality. Then it is easy to find integral subtours and add
the corresponding subtour elimination constraints to the ILP model. Iterating this pro-
cess until no more subtours are contained in the solution obviously solves the TSP to
optimality.
In this work we focus on the structure of subtour elimination constraints and how
to find a “good” set of subtour elimination constraints in reasonable time. Therefore,
we aim to identify the local structure of the vertices of a given TSP instance by run-
ning a clustering algorithm. Based on empirical observations and results from random
graph theory we further extend this clustering-based approach and develop a hierarch-
ical clustering method with a mechanism to identify subtour elimination constraints as
“relevant”, if they appear in consecutive iterations of the algorithm.
We mostly refrained from adding additional features which are highly likely to im-
prove the performance considerably, such as starting heuristics (cf. Section 4.4), lower
bounds or adding additional cuts. In the future it might be interesting to explore the
limits of performance one can reach with a purely integer linear programming approach
by adding these improvements. Clearly, we can not expect such a basic approach to com-
pete with the performance of Concorde [2], which has been developed over many years
and basically includes all theoretical and technical developments known so far. However,
it turns out that most of the standard benchmark instances with up to 400 vertices can
be solved in a few minutes by this purely integer strategy.
Finally, we briefly discussed some theoretical aspect for random Euclidean graphs
which could lead to polyhedral results in the expected case.

Acknowledgements
The research was funded by the Austrian Science Fund (FWF): P23829-N13.
We would like to thank the developers of the SCIP MIP-solver from the Konrad-
Zuse-Zentrum für Informationstechnik Berlin, especially Mr Gerald Gamrath, for their
valuable support.

References
[1] T. Achterberg, “Scip: Solving constraint integer programs,” Mathematical Program-
ming Computation, vol. 1, no. 1, pp. 1–41, Jul. 2009,

28
http://mpc.zib.de/index.php/MPC/article/view/4.

[2] D. L. Applegate, R. E. Bixby, V. Chvátal, and W. J. Cook, The Traveling Salesman

Problem: A Computational Study. Princeton University Press, 2006.

[3] J. Beardwood, J. H. Halton, and J. M. Hammersley, “The shortest path through

many points,” Mathematical Proceedings of the Cambridge Philosophical Society,
vol. 55, pp. 299–327, 1959.

[4] R. Bosch, “Connecting the dots: The ins and outs of TSP art,” in Bridges
Leeuwarden: Mathematics, Music, Art, Architecture, Culture. Winfield, Kansas:
Southwestern College, 2008, pp. 235–242.

[5] H. Crowder and M. W. Padberg, “Solving large-scale symmetric travelling salesman

problems to optimality,” Management Science, vol. 26, no. 5, pp. 495–509, 1980.

[6] G. Dantzig, R. Fulkerson, and S. Johnson, “Solution of a large-scale traveling-

salesman problem,” Operations Research, vol. 2, pp. 393–410, 1954.

[7] M. Grötschel and O. Holland, “Solving large-scale symmetric travelling salesman

problems to optimality,” Mathematical Programming, vol. 51, pp. 141–202, 1991.

[8] G. Gutin and A. Punnen, The Traveling Salesman Problem and Its Variations.
Springer, 2006.

[9] K. Helsgaun, “LKH – version 2.0.2,” website, 2008, available at

www.akira.ruc.dk/∼keld/research/LKH/LKH-2.0.2.tgz.

[10] A. K. Jain and R. C. Dubes, Algorithms for Clustering Data. Prentice Hall PTR,
1988.

[11] E. Lawler, J. Lenstra, A. Rinnooy Kan, and D. Shmoys, The Traveling Salesman
Problem: A Guided Tour of Combinatorial Optimization. J. Wiley, 1985.

[12] P. Miliotis, “Integer programming approaches to the travelling salesman problem,”

Mathematical Programming, vol. 10, pp. 367–378, 1976.

[13] D. Naddef and S. Thienel, “Efficient separation routines for the symmetric traveling
salesman problem ii: separating multi handle inequalities,” Mathematical Program-
ming, vol. 92, pp. 257–283, 2002.

[14] T. Öncan, İ. Kuban Altınel, and G. Laporte, “A comparative analysis of several
asymmetric traveling salesman problem formulations,” Computers & Operations
Research, vol. 36, no. 3, pp. 637–654, 2009.

[15] M. Padberg and G. Rinaldi, “An efficient algorithm for the minimum capacity cut
problem,” Mathematical Programming, vol. 47, pp. 19–36, 1990.

[16] M. Penrose, Random Geometric Graphs. Oxford University Press, 2003.

29
[17] G. Reinelt, “TSPLIB95,” website, 1995, available at
http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/.

[18] ——, The Traveling Salesman: Computational Solutions for TSP Applications.
Springer, 1994.

[19] W. T. Rhee, “A matching problem and subadditive Euclidean functionals,” The

Annals of Applied Probability, vol. 3, no. 3, pp. 794–801, 1993.

[20] A. Schrijver, Combinatorial Optimization: Polyhedra and Efficiency. Springer,

2003.

[21] R. Tošić and O. Bodroža, “On the number of hamiltonian cycles of P4 × Pn ,” Indian
Journal of Pure and Applied Mathematics, vol. 21, pp. 403–409, 1990.

[22] J. E. Yukich, Probability Theory of Classical Euclidean Optimization Problems.

Springer, 1998.

30
Appendix

Figure 5: Instance RE A 150. Euclidean Figure 6: Instance kroB150. Euclidean Figure 7: Instance u159. Euclidean dis-
distances between vertices. distances between vertices. tances between vertices.
31

Figure 8: Instance RE A 150: Main idea Figure 9: RE A 150: iteration 1. Figure 10: RE A 150: iteration 2.
of our approach – iteration 0.
Figure 11: RE A 150: iteration 3. Figure 12: RE A 150: iteration 4. Figure 13: RE A 150: iteration 5.
32

Figure 14: RE A 150: iteration 6. Figure 15: RE A 150: iteration 7. Figure 16: RE A 150: iteration 8.
Figure 17: RE A 150: iteration 9. Figure 18: RE A 150: iteration 10. Figure 19: RE A 150: iteration 11.

c′

27
33

24
21
18
15
12
9
6
3
n
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

Figure 20: Restricted clustering with c = n on random Euclidean graphs with minimum cluster size 3. The number of obtained
clusters c′ is plotted for every n. For every number of vertices n we created 100000 graphs.
t
100

2 without clustering (180.42 s)

c
0 5 10 15 20 25 30 35 40 45 50
best obtained time (16.54 s)
t.

20
34

10 without clustering (6.56 s)

without clustering (180.42 s)
c
0 5 10 15 20 25 30 35 40 45 50
best obtained time (3.8 s)

Figure 21: Computation time t in seconds depending on the number of clusters c for clustering (full line) and for restricted
clustering (dashed). Illustrative instances kroB150 (upper figure) and u159 (lower figure).
Mean number of iterations

12

n
0 35 70 105 140 175 210 245
Mean length

12

8
35

n
0 35 70 105 140 175 210 245

Figure 22: Mean number of iterations used by the BasicIntegerTSP (upper figure), mean length of an optimal TSP tour (lower
figure, dashed) and mean length of an optimal weighted 2-matching (lower figure, full line) in random Euclidean graphs. For
every number of vertices n we created 100 graphs.
Mean number of subtours

10
9
8
7
6
5
4
3
2
36

1
iteration ∗λ
0 1 2 3 4 5 6 7 8 9 10

Figure 23: Mean number of subtours during the BasicIntegerTSP in random Euclidean graphs for n = 60 sorted according
to the number of iterations (λ = 4/10 (full line), 5/10 (dashed), 6/10 (dotted), 7/10 (loosely dashed), 8/10 (loosely dotted)).
We created 100000 graphs.
instance BasicIntegerTSP HC | n HC | 4n/ log2 n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
ch150 13 7 74 150 2 435 9 5 223 14 6 129
kroA150 19 7 136 11 2 268 8 2 245 11 4 130
kroB150 179 8 148 72 2 301 34 4 315 21 4 168
pr152 16 13 184 5 3 214 5 3 205 9 4 174
u159 6 4 49 29 3 292 14 4 303 11 4 140
si175 52 10 183 99 6 494 40 8 415 44 7 263
brg180 44 4 103 54 2 185 102 18 359 24 2 27
rat195 347 6 274 241 3 491 272 4 419 267 5 322
d198 10986 10 301 483 7 894 1094 10 582 3986 9 326
kroA200 677 8 237 177 2 362 941 3 353 690 5 238
kroB200 31 5 121 37 3 292 23 3 269 31 4 164
gr202 39 11 77 2430 3 1216 61 8 330 60 8 217
tsp225 178 9 261 1113 3 981 138 6 551 151 6 341
37

pr226 5183 10 409 18 1 593 13 1 585 59 3 357

gr229 239 6 311 2984 4 1056 172 7 490 173 8 324
gil262 179 7 268 1052 2 807 169 3 564 217 4 368
a280 157 11 143 – – – 124 3 733 181 7 352
pr299 9263 9 413 4051 2 782 1998 5 745 1716 5 455
lin318 6885 8 357 457 2 756 274 5 660 275 5 355
rd400 2401 9 467 9329 5 1494 983 6 1018 1579 8 539
gr431 2239 9 453 – – – 4748 9 1734 4214 10 833
pcb442 2737 11 501 – – – 3830 16 1796 2277 15 888
u574 17354 6 423 – – – 18050 4 1290 8664 4 629
gr666 17711 8 789 – – – 23212 4 3104 18031 7 1408
rat783 30156 6 457 – – – – – – – – –
mean ratio 6.016088 0.890955 0.804727
RE A 150 23 8 100 103 5 363 36 5 238 28 6 162
RE B 150 13 7 78 99 1 424 13 3 255 14 4 146
instance BasicIntegerTSP HC | n HC | 4n/ log2 n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
RE C 150 9 5 70 21 1 235 7 3 195 7 3 98
RE D 150 8 6 60 50 2 274 7 4 197 7 4 112
RE E 150 9 7 55 76 1 339 22 4 274 17 5 149
mean RE 150 12.4 6.6 72.6 69.8 2 327 17 3.8 231.8 14.6 4.4 133.4
RE A 200 72 7 163 560 3 613 77 4 376 157 6 304
RE B 200 125 7 148 452 2 582 76 6 348 205 5 250
RE C 200 84 8 178 107 1 373 35 4 341 43 4 220
RE D 200 29 5 102 173 2 425 40 3 336 44 5 190
RE E 200 65 9 139 411 2 561 29 4 301 21 2 151
mean RE 200 75 7.2 146 340.6 2 510.8 51.4 4.2 340.4 94 4.4 223
RE A 250 138 9 186 1154 4 923 158 6 540 163 7 334
RE B 250 642 7 263 689 3 599 198 5 497 533 5 359
RE C 250 156 6 219 1545 1 846 57 3 333 136 4 275
38

RE D 250 273 6 220 501 2 542 192 6 479 199 5 292

RE E 250 103 5 156 339 1 675 70 4 511 105 4 252
mean RE 250 262.4 6.6 208.8 845.6 2.2 717 135 4.8 472 227.2 5 302.4
RE A 300 866 6 295 3574 2 1142 575 5 648 460 4 357
RE B 300 1411 8 348 4297 3 1059 627 4 672 865 6 402
RE C 300 1071 8 339 2071 3 848 236 6 567 687 7 474
RE D 300 229 6 290 2419 4 962 320 5 544 339 5 416
RE E 300 577 7 272 1543 3 726 283 6 526 436 6 344
mean RE 300 830.8 7 308.8 2780.8 3 947.4 408.2 5.2 591.4 557.4 5.6 398.6
RE A 350 411 5 252 3186 2 904 332 3 657 286 4 377
RE B 350 1021 8 339 11818 2 1102 1234 7 691 985 5 463
RE C 350 248 6 207 1243 2 936 232 4 750 358 5 390
RE D 350 1718 9 412 4271 3 1087 529 3 691 957 5 428
RE E 350 556 5 261 4560 3 1208 485 4 695 323 4 408
mean RE 350 790.8 6.6 294.2 5015.6 2.4 1047.4 562.4 4.2 696.8 581.8 4.6 413.2
RE A 400 8456 8 454 82054 3 1328 10463 5 980 8245 5 594
instance BasicIntegerTSP HC | n HC | 4n/ log2 n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
RE B 400 88849 7 438 1043519 2 1661 57469 5 879 39759 6 589
RE C 400 780 6 312 53580 3 1755 779 6 858 875 4 450
RE D 400 2052 8 451 122357 2 1998 1546 4 796 1081 4 434
RE E 400 2847 7 332 222902 2 1409 2450 4 660 2151 4 515
mean RE 400 20596.8 7.2 397.4 304882.4 2.4 1630.2 14541.4 4.8 834.6 10422.2 4.6 516.4
RE A 450 2107 5 279 – – – 2947 3 872 3595 4 535
RE B 450 68338 8 413 – – – 57921 8 906 94135 6 575
RE C 450 46360 10 596 – – – 33505 6 1166 13425 7 723
RE D 450 1212 6 368 – – – 1718 3 902 1644 4 520
RE E 450 1539 8 391 – – – 2438 5 1098 1345 7 637
mean RE 450 23911.2 7.4 409.4 – – – 19705.8 5 988.8 22828.8 5.6 598
RE A 500 15330 6 436 – – – 7531 4 1000 10323 5 629
RE B 500 16883 6 352 – – – 33464 4 1558 79362 5 727
39

RE C 500 3724 5 428 – – – 1337 4 955 1437 5 585

RE D 500 322951 9 567 – – – 339694 7 1236 307921 6 743
RE E 500 243378 9 679 – – – 110212 4 1113 134563 9 889
mean RE 500 120453.2 7 492.4 – – – 98447.6 4.6 1172.4 106721.2 6 714.6
mean ratio 12.132531 0.898420 1.040018
mean ratio all 9.760849 0.895621 0.951784
Table 8: Results for BasicIntegerTSP and for different variants of the approach which uses the hierarchical clustering. Mean
ratios refer to the arithmetic means over ratios between the running times of the particular approaches and the running time
of the BasicIntegerTSP. “sec.” is the time in seconds, “#i.” the number of iterations and “#c.” the number of subtour
elimination constraints added to the ILP before starting the last iteration. The entries “–” by TSPLIB instances cannot be
computed with 16 GB RAM or would take more than 12 hours.

• BasicIntegerTSP

• HC | n – hierarchical clustering; the constraints cannot be dropped and the maximum size of a solved cluster is u = n
(i.e. in fact, there is no upper bound)
• HC | 4n/ log2 n – hierarchical clustering; the constraints cannot be dropped and the maximum size of a solved cluster
is u = 4 logn n
2

• HCD | 4n/ log2 n – hierarchical clustering; the constraints can be dropped and the maximum size of a solved cluster
is u = 4 logn n
2
40
instance BasicIntegerTSP C | ⌊n/5⌋ RC3 | ⌊n/5⌋ RC3 | n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
ch150 13 7 74 12 6 114 9 6 109 16 7 117 14 6 129
kroA150 19 7 136 25 6 187 43 6 166 33 5 185 11 4 130
kroB150 179 8 148 53 4 219 138 5 215 44 5 202 21 4 168
pr152 16 13 184 17 12 181 17 11 204 18 12 181 9 4 174
u159 6 4 49 6 4 149 5 5 151 3 3 70 11 4 140
si175 52 10 183 31 10 213 55 13 250 35 9 196 44 7 263
brg180 44 4 103 17 3 81 19 8 102 121 11 316 24 2 27
rat195 347 6 274 246 4 268 275 6 315 114 6 257 267 5 322
d198 10986 10 301 4253 11 315 – – – 4762 9 321 3986 9 326
kroA200 677 8 237 332 6 214 350 5 190 287 4 171 690 5 238
kroB200 31 5 121 29 5 148 21 5 147 32 4 123 31 4 164
gr202 39 11 77 50 8 233 36 6 174 25 6 143 60 8 217
tsp225 178 9 261 100 9 223 84 10 235 100 8 300 151 6 341
41

pr226 5183 10 409 3614 6 363 36744 5 403 12944 9 415 59 3 357
gr229 239 6 311 335 6 289 152 6 256 311 7 341 173 8 324
gil262 179 7 268 250 8 250 133 7 268 152 6 274 217 4 368
a280 157 11 143 61 4 299 196 11 350 117 9 221 181 7 352
pr299 9263 9 413 6376 7 387 7410 6 416 16059 6 414 1716 5 455
lin318 6885 8 357 537 7 331 386 6 364 1560 6 391 275 5 355
rd400 2401 9 467 1212 7 420 1827 7 438 1522 8 398 1579 8 539
gr431 2239 9 453 3098 9 626 3384 9 647 2496 10 704 4214 10 833
pcb442 2737 11 501 3868 16 770 1815 17 567 2626 16 594 2277 15 888
u574 17354 6 423 11702 4 498 35204 5 580 13722 5 572 8664 4 629
gr666 17711 8 789 11756 7 919 14223 7 1001 13573 7 1002 18031 7 1408
rat783 30156 6 457 184381 5 701 37805 5 735 38630 6 779 – – –
mean ratio 1.014009 1.170299 0.983280 0.804727
RE A 150 23 8 100 18 5 142 26 6 141 36 7 155 28 6 162
RE B 150 13 7 78 8 4 117 13 5 129 7 4 86 14 4 146
instance BasicIntegerTSP C | ⌊n/5⌋ RC3 | ⌊n/5⌋ RC3 | n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
RE C 150 9 5 70 6 4 63 8 4 111 9 5 89 7 3 98
RE D 150 8 6 60 7 4 100 8 5 97 6 4 78 7 4 112
RE E 150 9 7 55 9 6 103 8 4 103 10 4 114 17 5 149
mean RE 150 12.4 6.6 72.6 9.6 4.6 105 12.6 4.8 116.2 13.6 4.8 104.4 14.6 4.4 133.4
mean RE 200 75 7.2 146 55.6 5.6 207.4 44.6 5.4 185.6 56.6 5.6 184.2 94 4.4 223
RE A 250 138 9 186 160 8 258 338 9 287 119 7 242 163 7 334
RE B 250 642 7 263 306 6 295 542 5 313 366 6 259 533 5 359
RE C 250 156 6 219 104 4 175 110 6 229 135 5 211 136 4 275
RE D 250 273 6 220 186 6 262 377 6 316 403 7 293 199 5 292
RE E 250 103 5 156 66 5 207 68 4 271 110 5 239 105 4 252
mean RE 250 262.4 6.6 208.8 164.4 5.8 239.4 287 6 283.2 226.6 6 248.8 227.2 5 302.4
RE A 300 866 6 295 1233 7 343 576 6 324 467 5 311 460 4 357
RE B 300 1411 8 348 1139 6 391 1100 7 431 1146 7 372 865 6 402
42

RE C 300 1071 8 339 608 6 392 331 6 314 458 6 312 687 7 474
RE D 300 229 6 290 276 6 321 268 7 307 396 7 374 339 5 416
RE E 300 577 7 272 353 7 322 353 5 320 464 6 334 436 6 344
mean RE 300 830.8 7 308.8 721.8 6.4 353.8 525.6 6.2 339.2 586.2 6.2 340.6 557.4 5.6 398.6
RE A 350 411 5 252 695 5 275 513 5 277 375 4 268 286 4 377
RE B 350 1021 8 339 793 7 363 1027 8 362 900 7 353 985 5 463
RE C 350 248 6 207 196 5 232 326 7 280 296 5 310 358 5 390
RE D 350 1718 9 412 749 8 385 1047 5 381 781 6 428 957 5 428
RE E 350 556 5 261 471 6 356 364 4 368 352 4 339 323 4 408
mean RE 350 790.8 6.6 294.2 580.8 6.2 322.2 655.4 5.8 333.6 540.8 5.2 339.6 581.8 4.6 413.2
RE A 400 8456 8 454 16648 6 471 24941 7 489 28803 7 516 8245 5 594
RE B 400 88849 7 438 72010 7 496 77325 7 503 59499 7 497 39759 6 589
RE C 400 780 6 312 1198 6 430 831 5 406 1095 7 453 875 4 450
RE D 400 2052 8 451 1639 5 436 591 5 454 1595 6 452 1081 4 434
RE E 400 2847 7 332 1602 6 434 3724 5 390 1608 6 408 2151 4 515
mean RE 400 20596.8 7.2 397.4 18619.4 6 453.4 21482.4 5.8 448.4 18520 6.6 465.2 10422.2 4.6 516.4
instance BasicIntegerTSP C | ⌊n/5⌋ RC3 | ⌊n/5⌋ RC3 | n HCD | 4n/ log2 n
sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c. sec. #i. #c.
RE A 450 2107 5 279 2921 4 333 2915 5 456 1385 4 383 3595 4 535
RE B 450 68338 8 413 15587 7 439 12828 5 442 24941 8 494 94135 6 575
RE C 450 46360 10 596 38930 9 697 35388 6 632 22898 7 647 13425 7 723
RE D 450 1212 6 368 948 6 460 2175 6 429 2120 7 388 1644 4 520
RE E 450 1539 8 391 2210 8 480 1786 7 434 1901 7 432 1345 7 637
mean RE 450 23911.2 7.4 409.4 12119.2 6.8 481.8 11018.4 5.8 478.6 10649 6.6 468.8 22828.8 5.6 598
RE A 500 15330 6 436 10907 6 576 14786 5 543 6118 6 531 10323 5 629
RE B 500 16883 6 352 12299 5 453 19681 4 483 186708 5 535 79362 5 727
RE C 500 3724 5 428 3063 6 519 2643 4 471 2339 5 440 1437 5 585
RE D 500 322951 9 567 514403 8 701 231961 6 618 314232 9 684 307921 6 743
RE E 500 243378 9 679 167194 8 718 125303 9 685 82051 9 671 134563 9 889
mean RE 500 120453.2 7 492.4 141573.2 6.6 593.4 78874.8 5.6 560 118289.6 6.8 572.2 106721.2 6 714.6
mean ratio 0.898743 0.964008 1.180427 1.040018
43

mean ratio all 0.943076 1.041367 1.104601 0.951784

Table 9: Comparison between different variants of our approach. Mean ratios refer to the arithmetic means over ratios between
the running times of the particular approaches and the running time of the BasicIntegerTSP. “sec.” is the time in seconds,
“#i.” the number of iterations and “#c.” the number of subtour elimination constraints added to the ILP before starting
the last iteration. The entries “–” by TSPLIB instances cannot be computed with 16 GB RAM.

• BasicIntegerTSP

• C | ⌊n/5⌋ – clustering for c = ⌊ n5 ⌋

• RC3 | ⌊n/5⌋ – restricted clustering for c = ⌊ n5 ⌋; the minimum size of a cluster is 3

• RC3 | n – restricted clustering for c = n; the minimum size of a cluster is 3

• HCD | 4n/ log2 n – hierarchical clustering; the constraints can be dropped and the maximum size of a solved cluster
is u = 4 logn n
2