Stochastic Runtime Analysis of a Cross-Entropy
Algorithm for Traveling Salesman Problems
Zijun Wua,1,∗, Rolf H. Möhringb,2,∗∗, Jianhui Laic,

arXiv:1612.06962v2 [] 11 Oct 2017

a

Beijing Institute for Scientific and Engineering Computing (BISEC), Pingle Yuan 100,
Beijing, P. R. China
b
Beijing Institute for Scientific and Engineering Computing (BISEC), Pingle Yuan 100,
Beijing, P. R. China
c
College of Metropolitan Transportation, Beijing University of Technology, Pingle Yuan
100, Beijing, P. R. China

Abstract
This article analyzes the stochastic runtime of a Cross-Entropy Algorithm
mimicking an Max-Min Ant System with iteration-best reinforcement. It
investigates the impact of magnitude of the sample size on the runtime to
find optimal solutions for TSP instances.
For simple TSP instances that have a {1, n}-valued distance function and
a unique optimal solution, we show that sample size N ∈ ω(ln n) results in a
stochastically polynomial runtime, and N ∈ O(ln n) results in a stochastically
exponential runtime, where “stochastically” means with a probability of 1 −
n−ω(1) , and n represents number of cities. In particular, for N ∈ ω(ln n),
we prove a stochastic runtime of O(N · n6 ) with the vertex-based random
solution generation, and a stochastic runtime of O(N · n3 ln n) with the edgebased random solution generation. These runtimes are very close to the best
known expected runtime for variants of Max-Min Ant System with best-sofar reinforcement by choosing a small N ∈ ω(ln n). They are obtained for
∗

Principal corresponding author
Corresponding author
Corresponding author
Email addresses: zijunwu@bjut.edu.cn (Zijun Wu), rolf.moehring@tu-berlin.de
(Rolf H. Möhring), laijianhui@bjut.edu.cn (Jianhui Lai)
1
The author is also affiliated with the School of Applied Mathematics and Physics at
Beijing University of Technology
2
The author is a professor emeritus of mathematics at Berlin University of Technology
∗∗

Preprint submitted to Theoretical Computer Science

October 12, 2017

the stronger notion of stochastic runtime, and analyze the runtime in most
cases.
We also inspect more complex instances with n vertices positioned on
an m × m grid. When the n vertices span a convex polygon, we obtain
a stochastic runtime of O(n4 m5+ ) with the vertex-based random solution
generation, and a stochastic runtime of O(n3 m5+ ) for the edge-based random
solution generation. When there are k ∈ O(1) many vertices inside a convex
polygon spanned by the other n − k vertices, we obtain a stochastic runtime
of O(n4 m5+ + n6k−1 m ) with the vertex-based random solution generation,
and a stochastic runtime of O(n3 m5+ + n3k m ) with the edge-based random
solution generation. These runtimes are better than the expected runtime
for the so-called (µ+λ) EA reported in a recent article, and again obtained
for the stronger notion of stochastic runtime.
Keywords: probabilistic analysis of algorithms, stochastic runtime analysis
of evolutionary algorithms, Cross Entropy algorithm, Max-Min Ant
System, (µ+λ) EA.
1. Introduction
The Cross Entropy (CE) algorithm is a general-purpose evolutionary algorithm (EA) that has been applied successfully to many N P-hard combinatorial optimization problems, see e.g. the book [1] for an overview. It
was initially designed for rare event simulation by Rubinstein [2] in 1997,
and thereafter formulated as an optimization tool for both continuous and
discrete optimization (see [3]).
CE has much in common with the famous ant colony optimization (ACO,
see [4]) and the estimation of distribution algorithms (EDAs, see [5]). They
all belong to the so-called model-based search paradigm (MBS), see [6]. Instead of only manipulating solutions, which is very typical in traditional
heuristics like Genetic Algorithms [7] and Local Search [8] and others, MBS
algorithms attempt to optimize the solution reproducing mechanism. In each
iteration, they produce new solutions by sampling from a probabilistic distribution on the search space. The distribution is often called a model in the
literature (see e.g. [6] and [9]). This model evolves iteratively by incorporating information from some elite solutions occurring in the search history,
so as to asymptotically model the spread of optimal solutions in the search

2

space. See the recent Thesis [9] for more details on MBS algorithms and
their mathematical properties.
An important issue for MBS algorithms is to determine a suitable size
for the sampling in each iteration. A large sample size makes each iteration
unwieldy, however a small sample size may mislead the underlying search
due to the randomness in the sampling. Sample size reflects the iterative
complexity (computational complexity in each iteration). Whether a large
sample size is harmful depends on the required optimization time (i.e., the
total number of iterations required to reach an optimal solution). This article
aims to shed a light on this issue by theoretically analyzing the relation
between sample size and optimization time for a CE variant that includes
also some essential features of the famous Max-Min Ant System (MMAS
[10]). To this end, a thorough runtime analysis is needed.
The theoretical runtime analysis of EAs has gained rapidly growing interest in recent years, see e.g. [11], [12], [13], [14], [15], [16], [17], [18], [19],
[20], [21], [22], [23], [24], and [25]. In the analysis, an oracle-based view of
computation is adopted, i.e., the runtime of an algorithm is expressed as
the total number of solutions evaluated before reaching an optimal solution.
Since the presence of randomness, the runtime of an EA is often conveyed
in expectation or with high probability. Due to the famous No Free Lunch
Theorem [26], the analysis must be problem-specific. The first steps towards
this type of analysis were made for the so-called (1+1) EA [11] on some test
problems that use pseudo boolean functions as cost functions, e.g., OneMax
[15], LeadingOnes [27] and BinVar [11]. Recent research addresses problems of practical importance, such as the computing a minimum spanning
trees (MST) [28], matroid optimization [29], traveling salesman problem [30],
the shortest path problem [23], the maximum satisfiability problem [31] and
the max-cut problem [32].
Runtime analysis generally considers two cases: expected runtime analysis
and stochastic runtime analysis. Expected runtime is the average runtime of
an algorithm on a particular problem, see, e.g., the runtime results of (1 + 1)
EA reported in [11]. Expected runtime reflects the oracle-based average performance of an algorithm. A mature technique for expected runtime analysis
is the so-called drift analysis [12]. However, this technique requires that the
algorithm has a finite expected runtime for the underlying problem. By [33],
drift analysis is not applicable to the traditional CE [3].
An algorithm with a smaller expected runtime need not be more efficient,
see [34] for details. In contrast, stochastic runtime provides a better under3

standing of the performance of a (randomized) EA. Stochastic runtime is a
runtime result conveyed with an overwhelming probability guarantee (see,
e.g., the classic runtime result of 1-ANT in [15]), where an overwhelming
probability means a probability tending to 1 superpolynomially fast in the
problem size. It therefore reflects the efficiency of an algorithm for most cases
in the sense of uncertainty. This article is concerned with stochastic runtime
analysis, aiming to figure out the relation between stochastic runtime and
magnitude of the sample size.
Runtime analysis of CE algorithms has been initiated in [33], where Wu
and Kolonko proved a pioneering stochastic runtime result for the traditional
CE on the standard test problem LeadingOnes. As a continuation of the
study of [33], Wu et al [34] further investigated the stochastic runtime of
the traditional CE on another test problem OneMax. The runtime results
reported in [33] and [34] showed that sample size plays a crucial role in
efficiently finding an optimal solution. In particular, Wu et al [34] showed
that if the problem size n is moderately adapted to the sample size N , then
4
the stochastic runtime of the traditional CE on OneMax is O(n1.5+ 3  ) for
arbitrarily small  > 0 and a constant smoothing parameter ρ > 0, which
beats the best-known stochastic runtime O(n2 ) reported in [13] for the classic
1-ANT algorithm, although 1-ANT employs a much smaller sample size (i.e.,
sample size equals one). Moreover, by imposing upper and lower bounds on
the sampling probabilities as was done in MMAS [10], Wu et al [34] showed
further that the stochastic runtime of the resulting CE can be significantly
improved even in a very rugged search space.
The present article continues the stochastic runtime analysis of [34], but
now in combinatorial optimization with a study of CE on the traveling salesman problem (TSP). We emphasize the impact of the magnitude of N on
the stochastic runtime, put ρ = 1, and consider a CE variant resembling an
MMAS with iteration-best reinforcement under two different random solution generation mechanisms, namely, a vertex-based random solution generation and an edge-based random solution generation.
Stochastic runtime for MMAS with iteration-best reinforcement on simple problems like OneMax has been studied in [20] and [25]. In particular,
Neumann et al [20] showed that to obtain a stochastically polynomial runtime for OneMax, N/ρ ∈ Ω(ln n) is necessary. We shall not only extend
this to TSP for the case of ρ = 1, but also prove that N ∈ ω(ln n) is already
sufficient to guarantee a stochastically polynomial runtime for simple TSP
instances.
4

TSP is a famous N P-complete combinatorial optimization problem. It
concerns finding a shortest Hamiltonian cycle on a weighted complete graph.
Existing algorithms exactly solving TSP generally have a prohibitive complexity. For instance, the Held-Karp algorithm [35] solves the problem with
a complexity of O(n2 2n ). A well-known polynomial time approximation algorithm for metric TSP is the so-called Christofides algorithm [36], which
finds a solution with a cost at most 3/2 times the cost of optimal solutions.
As mentioned in [37], this is still the best known approximation algorithm
for the general metric TSP so far. For Euclidean TSP there exists a famous
polynomial-time approximation scheme (PTAS) by Arora, see [38]. To design a superior approximation algorithm, researchers in recent years tend to
study TSP instances with particular structures, see, e.g., [39].
Due to the prohibitive running time of exact algorithms, heuristics are
frequently employed in practice so as to efficiently compute an acceptable
solution for a TSP problem, e.g., the Lin-Kernighan (LK) algorithm [40].
As a popular heuristic, CE has also been applied in practice to solve TSP
instances, see [41] and [3]. The implementation there shows that CE can also
efficiently compute an acceptable solution.
In view of the high complexity of general TSP, we consider in our analysis
two classes of TSP instances with a particular structure. The first kind of
instances has been used in [19] and [42] for analyzing the expected runtime of
some MMAS variants with best-so-far reinforcement. These TSP instances
have polynomially many objective function values and a unique optimal solution. Moreover, on these TSP instances, solutions containing more edges
from the optimal solution have a smaller cost than those with fewer such
edges. For more details on these instances, see Section 5.
For these simple instances, we prove in Theorem 2 that with a probability

1 − e−Ω(n ) , the runtime is O(n6+ ) with the vertex-based random solution
generation, and O(n3+ ln n) with the edge-based random solution generation,
for any constant  ∈ (0, 1) and N ∈ Ω(n ). For the case of N ∈ ω(ln n), we
show that the runtimes (resp., O(n6 N and n3 (ln n)N ) are even smaller with
probability 1 − n−ω(1) , see Corollary 1. These results are very close to the
known expected runtime O(n6 + n lnρ n ) for (1 + 1) MMAA reported in [19],
and the expected runtime O(n3 ln n + nρln ) for MMAS∗Arb reported in [42]
(where ρ ∈ (0, 1) is an evaporation rate), if N ∈ ω(ln n) is suitably small.
But they give the stronger guarantee of achieving the optimal solution in the
respective runtime with an overwhelming probability. Moreover, we show

5

a stochastically exponential runtime for a suitable choice of N ∈ O(ln n),
see Theorem 3. This generalizes the finding in [20] for OneMax to TSP
instances. Therefore, N ∈ Ω(ln n) is necessary, and N ∈ ω(ln n) is sufficient
for a stochastically polynomial runtime for simple TSP instances.
We also inspect more complex instances with n vertices positioned on an
m×m grid, and the Euclidean distance as distance function. These instances
have been employed in [43] and [30] for analyzing the expected runtime of
(µ+λ) EA and randomized local search (RLS). When the n vertices span a
convex polygon without vertices in the interior of the polygon (so they are
the corners of that polygon), we prove a stochastic runtime of O(n4 m5+ )
for the vertex-based random solution generation, and a stochastic runtime
of O(n3 m5+ ) for the edge-based random solution generation, see Theorem
4 for details. Similarly, the  in the stochastic runtimes can be removed
by slightly decreasing the probability guarantee, see Corollary 2. When the
vertices span a convex polygon with k ∈ O(1) vertices in the interior, we
show a stochastic runtime of O(n4 m5+ + n6k−1 m ) with the vertex-based
random solution generation, and a stochastic runtime of O(n3 m5+ + n3k m )
with the edge-based random solution generation, see Theorem 5 for details.
These runtimes are better than the expected runtime for the so-called (µ+λ)
EA and RLS reported in the recent paper [30].
The remainder of this paper is arranged as follows. Section 2 defines the
traditional CE and related algorithms, Section 3 defines the traveling salesman problem and provides more details of the used CE variants, Section 4
shows some important facts on the two random solution generation methods,
and Section 5 reports the stochastic runtime results on the TSP instances.
A short conclusion and suggestions for future work are given in Section 6.
Notations for runtime
Our analysis employs some commonly used notations from complexity
theory. We use O(f (n)) to denote the class of functions which are bounded
from above by the function f (n), i.e., those functions g(n) with g(n) ≤ c·f (n)
for large enough n and some constant c ≥ 0 not depending on n. Similarly,
Ω(f (n)) is the class of functions that are bounded from below by f (n), i.e.,
for any g(n) ∈ Ω(f (n)) there exists a constant c > 0 not depending on n
such that g(n) ≥ c · f (n) for large enough n. Class Θ(f (n)) is the intersection of Ω(f (n)) and O(f (n)). Class o(f (n)) is the class of functions g(n)
with g(n)/f (n) → 0 as n → ∞, and class ω(f (n)) is the class of functions
6

g(n) with g(n)/f (n) → +∞ as n → ∞. Obviously, o(f (n)) ⊂ O(f (n)) and
ω(f (n)) ⊂ Ω(f (n)).
2. The general cross entropy algorithm and related algorithms
We now introduce the traditional CE algorithm. The CE variant we will
analyze inherits the framework of this traditional version. To compare our
results with those in the literature, we shall give also details about some
related algorithms.
2.1. The traditional cross entropy algorithm
Algorithm 1 lists the traditional CE that was proposed in [3], adapted to
an abstract notion of combinatorial optimization problems. The algorithm
assumes a combinatorial minimization problem (S, f ), where S is a finite
search space of “feasible” solutions and f is the cost function. Every feasible
solution s ∈ S is composed of elements from a fixed finite set A, the ground
set of the problem, i.e., we assume S ⊆ An for some integer n ∈ N. Furthermore there is a product distribution on the product space An that induces a
distribution on S ⊆ An . The distribution on An can usually be represented as
a vector (or matrix) of real-valued probabilities. The convex combination of
the two distributions in Step 6 of Algorithm 1 then corresponds to a convex
combination of the two vectors (or matrices).
Specific to the TSP, the ground set A can be the set of nodes or edges, n
is the number of nodes, and a feasible solution is a sequence of elements from
A that forms a Hamiltonian cycle. The product distribution for the TSP is
represented as an n×n matrix.
When we consider the set of nodes as our ground set A, each row i
of the matrix is a marginal distribution that specifies choice probabilities
for all nodes following the current node i. A random Hamiltonian cycle
is sequentially constructed from the product distribution by allowing only
nodes not yet visited as continuations in each step, see Algorithm 2 for more
details.
When we consider the set of edges as A, marginals of the product distribution will be represented by the same n×n matrix where the sum of the
(i, j)-th and (j, i)-th entries reflects the probability that the edge {i, j} occurs in a random solution. A random Hamiltonian cycle is still constructed
sequentially and only edges leading to a feasible solution are taken in each
step, see Algorithm 3 for details.
7

Algorithm 1 The general Cross-Entropy algorithm
Input:
an initial distribution Π0 on the solution space, a fixed smoothing parameter ρ ∈ (0, 1], a sample size N ∈ N+ , an elite size M ∈ N+ with M ≤ N
1:
2:
3:
4:
5:
6:
7:
8:

t = 0;
loop
(1)
(N )
independently generate N random solutions Xt , . . . , Xt with the
current distribution Πt ;
[1]
sort these N solutions in non-decreasing order as f (Xt ) ≤ · · · ≤
[N ]
f (Xt ) according to the cost function f ;
learn an empirical distribution Wt from the M best solutions
[1]
[M ]
X t , . . . , Xt ;
set Πt+1 = (1 − ρ)Πt + ρWt ;
t = t + 1;
end loop

Traditionally, CE sets a small elite ratio α ∈ (0, 1) and uses the best
bα · N c solutions in Step 5 to build the empirical distribution Wt . Here, we
use the elite size M instead. This does not intrinsically change the original
algorithm. Steps 3 and 5 depend on the detailed definition of the underlying
problem. We shall give details to them in Subsection 3.2.
Step 6 of Algorithm 1 plays a crucial role in the different theoretical analyses of the algorithm, see, e.g., [44], [33], [45], [9], [34]. The occurrence of
good solutions are probabilistically enforced by incorporating the new information Wt into Πt+1 . This idea, somehow, coincides with the reinforcement
learning in [46]. The smoothing parameter ρ reflects the relative importance
of the new information Wt in the next sampling. It balances global exploration and local exploitation to a certain degree. A larger ρ makes the
algorithm concentrate more on the particular area spanned by the elite solu[1]
[M ]
tions Xt , . . . , Xt , while a smaller ρ gives more opportunities to solutions
outside that area.
However, balancing global exploration and local exploitation through tuning ρ is ultimately limited. Wu and Kolonko [33] proved that the famous
“genetic drift” [47] phenomenon also happens in this algorithmic scheme,
i.e., the sampling (Step 3) eventually freezes at a single solution and that
solution needs not to be optimal. This means that the algorithm gradually
8

loses the power of global exploration.
As a compensation for global exploration, Wu et al [34] proved that a
moderately large sample size N might be helpful. The results there showed
that a moderately large N configured with a large ρ (e.g., ρ = 1) can make the
algorithm very efficient. Although a large N introduces a high computational
burden in each iteration, the total number of iterations required for getting
an optimal solution is considerably reduced.
Wu et al [34] also indicated another way to compensate the global exploration, i.e., imposing a lower bound πmin ∈ (0, 1) and an upper bound
πmax ∈ (0, 1) on the sampling distributions in each iteration. This idea is
originated from MMAS [10]. In each iteration t, after applying Step 6, the
entries of distribution Πt+1 that are out of the range [πmin , πmax ] are reset to
that range by assigning to them the nearest bounds, see (6) in Section 3 for
more details. Wu et al [34] have proved that this can make CE more efficient
even in the case of a rugged search space.
To follow these theoretical suggestions made in [34], we shall in our
stochastic runtime analysis use a CE that modifies the traditional CE (Algorithm 1) accordingly. We shall see that these modifications make the CE
very efficient for the considered TSP instances.
2.2. Related evolutionary algorithms
Related evolutionary algorithms for TSP whose runtime has been extensively studied are RLS [28], (µ + λ) EA [30], and those theoretical abstractions of MMAS [10] including MMAS∗bs [17], (1+1) MMAA [19]. We now
give algorithmic details of them. In order to facilitate the comparison, their
runtimes for TSP instances will be discussed in Section 5.
(µ + λ) EA is an extension of the famous (1 + 1) EA [11]. (µ + λ) EA
randomly chooses µ solutions as the initial population. In each iteration,
(µ+λ) EA randomly chooses λ parents from current population, then produces
λ children by applying randomized mutation to each of the selected parents,
and forms the next population by taking the best µ solutions from these
µ + λ solutions in the end of current iteration. The expected runtime of
(µ + λ) EA on TSP instances is studied in [30], where Sutton et al uses a
Poisson distribution to determine the number of randomized mutations (2opt move or jump operation) should be taken by a selected parent in each
iteration.
RLS is a local search technique [48]. It employs a randomized neighborhood. In each iteration, it randomly chooses a number of components of the
9

best solution found so far and then changes these components. The expected
runtime of RLS for TSP instances is also studied in [30], where the neighborhood is taken to be a k-exchange neighborhood with k randomly determined
by a Poisson distribution.
(1+1) MMAA is a simplified version of the famous MMAS [10], where
the sample size is set to 1 and pheromones are updated only with the best
solution found so far (best-so-far reinforcement) in each iteration. In each
iteration of (1+1) MMAA, the ant which constructed the best solution found
so far deposits an amount πmax of pheromones on the traversed edges, and an
amount πmin of pheromones on the non-traversed edges, and the pheromones
are updated by linearly combining the old and these newly added pheromones
as in Algorithm 1. The expected runtime of (1+1) MMAA on simple TSP
instances is studied in [19]. The expected runtime of its variant MMAS∗Arb
on simple TSP instances is studied in [42].
3. The traveling salesman problem and details of the CE variant
Now, we formally define TSP, and give more details of the CE variant we
will analyze.
3.1. The traveling salesman problem
We consider an undirected
graph G = (V, E) with vertex set V = {1, . . . , n}

and edge set E = {i, j} | i ∈ V, j ∈ V, i 6= j . A Hamiltonian cycle is a
sequence {{il , il+1 } | l = 1, . . . , n} of edges such that
a) i1 = in+1 ;
b) (i1 , . . . , in ) is a permutation of {1, 2, . . . , n}.
This definition actually considers E as the ground set A. As mentioned above,
we can also put A = V and represent Hamiltonian cycles in a more compact way as permutations of V. Note that a Hamiltonian cycle corresponds
to n different permutations, whereas a permutation corresponds to a unique
Hamiltonian cycle. However, the two representations are intrinsically the
same. We shall use them interchangeably in the sequel. To facilitate our
discussion, we shall refer to a Hamiltonian cycle by just referring to one of
the n corresponding permutations, and denote by S the set of all possible
permutations. We employ the convention that two permutations are said

10

to be same iff they form the same underlying Hamiltonian cycle. The notation {k, l} ∈ s shall mean that the edge {k, l} belongs to the underlying
Hamiltonian cycle of the solution (permutation) s.
Once a distance function d : E 7→ R+ is given, the (total traveling) cost
f (s) of a feasible solution s = (i1 , i2 , . . . , in ) ∈ S is then calculated by
f (s) :=

n−1
X

d(ij , ij+1 ) + d(in , i1 ).

(1)

j=1

We denote by S ∗ ⊆ S the set of feasible solutions (Hamiltonian cycles) that
minimize the cost (1).
3.2. Details of the CE variant
The CE variant we consider in the analysis completely inherits the structure of Algorithm 1, and additionally employs a component from MMAS.
We now formalize the sampling distribution, and define Steps 3 and 5 in
more detail. As mentioned, we represent a sampling distribution (a product
distribution on An ) for the TSP by a matrix Π = (πi,j )n×n , such that
Pn
a)
j=1 πi,j = 1, for all i = 1, . . . , n,
b) πi,i = 0 for all i = 1, . . . , n,
c) πi,j = πj,i for each edge {i, j} ∈ E.
For each edge {i, j} ∈ E, πi,j reflects the probability that a Hamilton cycle
continues with vertex j when it is in vertex i. In the sequel, we write the
t
sampling distribution Πt in iteration t as (πi,j
)n×n , where the superscript t of
t
0
πi,j indicates the iteration. The initial distribution Π0 = (πi,j
)n×n is, without
1
0
0
for
loss of generality, set to be the uniform distribution, i.e., πi,j = πj,i
= n−1
all edges {i, j} ∈ E.
We shall consider two random solution generation methods, a vertex-based
random solution generation and an edge-based random solution generation.
Algorithm 2 lists the vertex-based random solution generation method. This
method uses V as the ground set A. A product distribution of An is therefore
represented as a matrix Π = (πi,j )n×n satisfying a)-c) above, i.e., each row
of Π represents a sampling distribution on A = V. Directly sampling from
Π may produce infeasible solutions from An − S. To avoid that, Algorithm
2 starts with a randomly fixed initial node, and then sequentially extends
11

a partial solution with an unvisited vertex until a complete permutation is
obtained. This method is efficient and rather popular in practice, see, e.g.,
[41] and [4]. Here, “s + (v)” means that appends a vertex v to the end of a
partial solution s.
Algorithm 2 Vertex-based random solution generation
Input:
a distribution Π = (πi,j )n×n
Output:
a permutation of 1, 2, . . . , n
1: s = ∅, and Vunivisted = V ;
2: randomly select v from V, s = s + (v), and Vunvisited = Vunvisited − {v};
3: while (|Vunvisited 6= ∅|) do
4:
select a random vertex v 0 from Vunvisited with a probability
P[v 0 | s] = P

πv,v0

k∈Vunvisited

5:
6:
7:
8:

πv,k

;

(2)

set s = s+(v 0 ), Vunvisited = Vunvisited − {v 0 };
v = v0;
end while
return s;

The edge-based random solution generation is listed in Algorithm 3. The
idea is from [42]. This method considers edge set E as the ground set A. A
feasible solution is then a sequence of edges that form a Hamiltonian cycle, i.e.
S ⊆ E n . To unify the notation of feasible solutions, Algorithm 3 translates
its outcomes into permutations. As the actual ground set is E, a product
matrix such that each row is a marginal specifying
distribution is an n× n(n−1)
2
a sampling distribution on E. Algorithm 3 only considers those with identical
marginals, a product distribution can be therefore fully characterized by
one of its marginals and is therefore again represented by an n × n matrix
Π = (πi,j )n×n as above. PAn edge
P {i, j} ∈ E is then sampled from Π with
probability (πi,j + πj,i )/ nk=1 nl=1 πk,l = 2πi,j /n since each row of Π sums
up to 1. A random sequence ∈ E n is generated by independently sampling
from Π n times. To avoid infeasible solutions, Algorithm 3 considers in
every sampling only edges that are admissible by the edges selected before.
Given a set B of edges such that the subgraph (V, B) does neither contain a
12

cycle nor a vertex of degree ≥ 3, an edge e0 ∈ E is said to be admissible by
B if and only if the subgraph (V, B ∪ {e0 }) still does neither contain a cycle
nor a vertex of degree ≥ 3. We denote by Badmissible the set of edges ∈
/ B that
are admissible by B.
Algorithm 3 Edge-based random solution generation
Input:
a distribution Π = (πi,j )n×n
Output:
a permutation of 1, 2, . . . , n
1: B = ∅, Badmissible = E;
2: while (|B| ≤ n − 1) do
3:
select an edge {i, j} from Badmissible with a probability
P[e | s] = P
4:
5:
6:
7:
8:

πi,j + πj,i
;
{k,l}∈Badmissible πk,l + πl,k

(3)

set B = B∪{{i, j}};
update Badmissible ;
end while
let s = (1, i2 , i3 , . . . , in ) with {1, i2 }, {ij , ij+1 } ∈ B for j = 2, . . . , n−1;
return s;
(1)

(N )

The N random solutions Xt , . . . , Xt in iteration t are then generated
by N runs of Algorithm 2 or Algorithm 3 with the current distribution Πt =
t
t
)n×n is then calculated from
)n×n . The empirical distribution Wt = (wi,j
(πi,j
the M elite solutions by setting
PM
({i, j})
k=1 1{e0 ∈E | e0 ∈X[k]
t
t }
wi,j
=
,
(4)
M
where 1A (·) is the indicator function of set A = {e0 ∈ E | e0 ∈ Xt } for each
t+1
{i, j} ∈ E. The next distribution Πt+1 = (πi,j
)n×n is therefore obtained as
[k]

t+1
t
t
πi,j
= (1 − ρ)πi,j
+ ρwi,j

(5)

for each {i, j} ∈ E.
We continue with the suggestions made in [34]. In the CE variant, we
shall use a moderately large N and a large ρ = 1. To fully use the best elite
13

solutions, we take M = 1. To prevent premature convergence (i.e., a possible
stagnation at a non-optimal solution), we employ a feature from MMAS
[10], called max-min calibration, in the construction of Πt+1 ,. We choose
a lower bound πmin ∈ (0, 1) and an upper bound πmax ∈ (0, 1), and, after
applying (5), adjust Πt+1 by

t+1

πmin if πi,j < πmin ,
t+1
t+1
t+1
πi,j
= πi,j
(6)
if πi,j
∈ [πmin , πmax ],


t+1
πmax if πi,j
> πmax ,
for any edge {i, j} ∈ E. Note that the max-min calibration is the only step
that does not occur in the general CE (i.e., Algorithm 1).
This setting turns CE into an MMAS with iteration-best reinforcement,
[1]
i.e., only the iteration-best solution Xt is allowed to change the ‘pheromones’
Πt . Stützle and Hoos [10] indicated in an empirical study that the practical performance of iteration-best reinforcement is comparable to best-so-far
reinforcement for TSP instances. Thus, it should also be worthwhile to compare the theoretical runtime of iteration-best reinforcement with the known
expected runtimes of best-so-far reinforcement for TSP instances presented
in, e.g., [19] and [42].
4. Properties of the random solution generation methods
Before we start with our runtime analysis, we shall discuss some relevant
properties of the two random solution generation methods, which concern the
probability of producing a k-exchange move of the iteration-best solution in
the next sampling.
Formally, a k-exchange move on a Hamiltonian cycle is an operation that
removes k edges from the cycle and adds k new edges to obtain again a cycle.
A k-opt move is a k-exchange move reducing the total travel cost. Figure
1a shows an example of a 2-exchange move, in which edges {i, j}, {k, l} are
removed, and edges {i, l}, {k, j} are added.Figure 1b shows an example of a
3-exchange move.
In our analysis, we shall consider only iteration-best reinforcement with
ρ = 1 and the max-min calibration (6). The empirical distribution Wt =
t
(wi,j
)n×n for each iteration t ∈ N in this particular case therefore satisfies
(
[1]
min{1, πmax } = πmax if edge {i, j} ∈ Xt ,
t+1
t+1
πi,j
= πj,i
=
(7)
max{0, πmin } = πmin otherwise,
14

i

k

i

l

w

v

k

j

l

j
(a) A 2-exchange move

(b) A 3-exchange move

Figure 1: Examples for edge exchange moves

for every edge {i, j} ∈ E and iteration t ∈ N. Furthermore, Πt+1 = Wt .
[1]
Since Πt+1 is biased towards the iteration-best solution Xt , k-exchanges
[1]
of Xt with a large k are unlikely to happen among the N draws from Πt+1
by either of the two generation methods. Thus, an optimal solution is more
likely to be reached by a sequence of repeatedly k-exchange moves with small
k from iteration-best solutions. Therefore, it is necessary to estimate the
[1]
probabilities of producing a k-exchange of Xt in the two generation methods,
especially for the case of small k.
4.1. Probabilities of producing k-exchanges in the vertex-based random solution generation
The probability of producing k-exchanges with k = 2, 3 in the vertexbased random solution generation has been studied in Zhou [19]. With πmin =
1
and πmax = 1 − n1 , Zhou [19] proved for (1 + 1) MMAA that with a
n2
probability of Ω(1/n5 ), Algorithm 2 produces a random solution having more
edges from s∗ than x∗t (the best solution found so far) provided that x∗t is not
optimal. Zhou [19] actually showed that if x∗t 6= s∗ , then there exists either
a 2-opt move or a 3-opt move for x∗t , and Algorithm 2 produces an arbitrary
2-exchange of x∗t with a probability of Ω(1/n3 ), and an arbitrary 3-exchange
of x∗t with a probability of Ω(1/n5 ).
1
Although we use πmin = n(n−2)
and consider iteration-best reinforcement,
a similar result holds in our case. Claim 1 below gives a lower bound on the
probability of producing a k-exchange move of the iteration-best solution in
the next round with the vertex-based random solution generation.

15

[1]

Claim 1. Let M = 1, ρ = 1, and consider a k-exchange move of Xt for some
integer k = 2, 3, . . . , n. Then, Algorithm 2 produces the given k-exchange
move with a probability Ω(1/n2k−1 ) in every of the N draws in iteration t + 1.
Proof. Recall that in Algorithm 2, the probability (2) to select a continuing
t
t
) for
(or, equivalently, πj,i
edge {i, j} is always bounded from below by πi,j
each iteration t ∈ N, since each row of Πt sums up to 1. Given a k-exchange
[1]
move of Xt , one possibility to generate it from Πt+1 = Wt by Algorithm 2
is that one of the new edges is added in the last step. This happens with a
n−k
 1 k−1
1
· 1 − n1
≥ e·n2k−1
probability at least n1 · n(n−2)
, where e ≈ 2.71828
1
is Euler’s number, n represents the probability to select the starting vertex,
1
is the common lower bound of the probability to select the remaining
n(n−2)
k − 1 new edges, and 1 − n1 is the common lower bound of the probability to
[1]
select one of the remaining n − k edges from Xt .
[1]

Because of Claim 1, every 2-exchange of Xt is produced from Πt+1 by
Algorithm 2 with a probability Ω(1/n3 ), and every 3-exchange is produced
by Algorithm 2 with a probability Ω(1/n5 ). Note that for any k = 2, 3, . . . ,
[1]
if a k-opt move of Xt occurs among the N draws in the next sampling,
[1]
[1]
then f (Xt+1 ) < f (Xt ) must hold. Thus, if we take a moderately large
sample size, say N = Θ(n5+ ) for some  > 0, with a probability 1 − (1 −

[1]
[1]
Ω(n−5 ))Ω(5+) = 1 − e−Ω(n ) , f (Xt+1 ) < f (Xt ) will hold, provided that there
[1]
still exists a 2-opt or 3-opt move from Xt .
Claim 2. Suppose that M = 1, ρ = 1. Then, for iteration t + 1, the probability
[1]
that Algorithm 2 produces a solution with a cost not larger than Xt in one
application is in Ω(1).
[1]

Proof. Observe that the probability that Xt is reproduced in one application
of Algorithm 2 is larger than (1 − 1/n)n−1 ∈ Ω(1), which implies that the
[1]
cost of the generated random solution is not larger than f (Xt ).
[1]

Note that if Xt is reproduced at least once among the N draws in the
[1]
next sampling, then f (Xt+1 )[1] ≤ f (Xt ). Thus, if the sample size N ∈
[1]
[1]
Ω(ln n), then f (Xt+1 ) ≤ f (Xt ) with a probability 1 − (1 − Ω(1))N = 1 −
[1]
[1]
O(1/n). Particularly, when N ∈ Ω(n ) for some  > 0, f (Xt+1 ) ≤ f (Xt )

with an overwhelming probability 1 − e−Ω(n ) .
16

4.2. Probabilities of producing k-exchanges in the edge-based random solution
generation
The behavior of the edge-based random solution generation is comprehensively studied in [42]. Kötzing et al [42] proved for MMAS∗Arb and a constant
k ∈ O(1) that, with a probability of Ω(1), Algorithm 3 produces a random
solution that is obtained by a k-exchange move from the best solution found
so far.
t
t
t
t
Recall that in each iteration t, either πi,j
= πj,i
= πmin or πi,j
= πj,i
= πmax
for any edge {i, j} ∈ E. For convenience, we will call an edge {i, j} ∈ E with
t
t
= πmax a high edge, and otherwise a low edge. Kötzing et al [42]
= πj,i
πi,j
showed the probability of the event that Algorithm 3 chooses
√ a high edge in
an arbitrary fixed step conditioned on the event that l ≤ n low edges have
been chosen in some l steps before this step is 1 − O(1/n). Our setting is
1
only slightly different with from theirs, i.e., we use πmin = n(n−2)
but they
1
put πmin = n(n−1) . Thus, the result should also hold here. Claim 3 below
formally asserts this, readers may also refer to [42] for a similar proof.
Claim 3. Assume M = 1, ρ = 1. Then, the probability of choosing√a high
edge at any fixed step in Algorithm 3 is at least 1 − 12/n if at most n low
edges have been chosen before that step and there exists at least one high
admissible edge to be added.
Proof.
√ We now fix a step n−m for some m = 0, 1, . . . , n−1, and assume that
l ≤ n low edges have been chosen before this step. Obviously, we still need
to add m + 1 ≥ 1 edges to obtain a complete solution. We now estimate the
numbers of admissible high and low edges in this step. Note that every of
the l low edges blocks at most 3 of the m+l remaining high edges (at most
two which are incident to the end points of the low edge, and at most one
that may introduce a cycle). So at least m+l − 3l = m−2l ≥ m−3l high
edges are available for adding in this step. Of course, it may happen that
there is no admissible high edges in this step. However, we are not interested
in such a case. We consider only the case that there exists at least one
admissible high edge in this step, i.e. the number of admissible high edges
in this step is at least max{1, m − 3l}. Note also that the n−m edges added
before partition the subgraph of G = (V, E) with vertices V and edges from
the partial solution constructed so far into exactly m connected components
(here, we see an isolated vertex also as a connected component). For any two
of the components, there are at most 4 admissible edges connecting them.
17

 
Therefore,
there are at most min{4 m2 , n2 } admissible low edges. Observing
√
l ≤ n, the probability of choosing a high edge in this step is bounded from
below by
(
 
√
2m2
3
min{4 m2 , n2 } πmin
1− (m−3l)n(n−2)
≥ 1− (n−2)
if m > 3 n,
≥
(8)
1−
√
12
max{1, m − 3l} πmax
if m ≤ 3 n,
1− n−2
where the first inequality is obtained by observing that
   
 
m
n
n
2
min{4
,
} ≤ min{2m ,
} ≤ 2m2 ,
2
2
2
τmax = 1 − 1/n, πmin =

1
,
n(n−2)

and

2m2
m−3l

≤

2

√
1
− 3 2n
m
m

≤ 3n.

With Claim 3, we can show that, for any t ∈ N and any fixed k ∈ O(1),
[1]
the probability of the event that a k-exchange of Xt is produced by one
application of Algorithm 3 is Ω(1), see Claim 4. Here, we shall use a different
proof from the one presented by Kötzing et al [42], which appears to us as
problematic.
Claim 4. Let M = 1, ρ = 1. For any k ∈ O(1), with probability Ω(1), the
[1]
random solution produced by Algorithm 3 is a k-exchange of Xt .
Proof. Let k ∈ O(1) be arbitrarily fixed, and M be the set of all k-element
subsets of {1, 2, . . . , n/2} (where we assume without loss of generality that
n is even). Obviously, |M| ∈ Θ(nk ) since k ∈ O(1). Let M ∈ M be an arbitrarily fixed k-element subset. The probability of the event that Algorithm 3
selects k new edges (low edges) at steps i ∈ M and n − k edges (high edges)
[1]
from Xt at other steps, is bounded from below by

n−i+1
− (n − i + k + 1))πmin
1
1 n−k Y ( 2
≥ Θ( k ), (9)
1 − O( )
n
n(n − 1)πmin + (n − i + k + 1)πmax
n
i∈M
where 1 − O(1/n) is a lower bound for the probability of selecting an edge
[1]
from Xt . In each step i ∈ M, the edges chosen before partition the graph
into n−i+1 connected components, and for any two of the components there
exists at least 2 edges connecting
them without introducing a cycle. Hence,

n−i+1
there are at least
admissible edges in each step i ∈ M. Notice also
2
that the number of admissible high edges in this case is at most n − i + k + 1
18

(n − i+k+1 is the maximal number of high edges that have not been chosen
((n−i+1
)−(n−i+k+1))πmin
2
of (9) is just the lower
before). Therefore, each factor n(n−1)π
min +(n−i+k+1)πmax
bound of the probability for choosing an admissible edge not belonging to
[1]
Xt in a step i ∈ M.
As a result, the probability of the random event that Algorithm 3 pro[1]
duces a k-exchange of Xt with k ∈ O(1) in any of the N independent draws
in iteration t+1 is bounded from below by |M|·Θ( n1k ) = Θ(nk )·Θ( n1k ) ∈ Ω(1),
since new edges can also be added in steps l ≥ n/2.
Notice that in the edge-based random solution generation, for any k =
[1]
2, 3, . . . , n, any two k-exchanges of Xt are generated with the same probability, since the generation does not require adding the edges in a particular
order. Therefore, by Claim 4, for any k ∈ O(1), any specified k-exchange of
[1]
[1]
Xt will be produced with a probability Θ(1/nk ). Since reproducing Xt can
[1]
be seen as a 0-exchange of Xt , we can thus derive the following conclusion.
Claim 5. Let M = 1, ρ = 1. With probability Ω(1), the random solution
[1]
generated by Algorithm 3 has a cost not larger than that of Xt .
Claim 6 shows that it is unlikely that the random solution generated by
[1]
Algorithm 3 is “very” different from the last iteration-best solution Xt . This
will be fundamental for deriving the runtime lower bound.
Claim 6. Let M = 1, ρ = 1. For any δ ∈ (0, 1], with an overwhelming
min{δ,1/4}/2 )
probability 1 − e−ω(n
, the random solution generated by Algorithm 3
[1]
is a k-exchange move from Xt for some k < nδ .
Proof. Let δ ∈ (0, 1] be arbitrarily fixed, and put γ = min{δ, 1/4}. To prove
the claim, we just need to show that with an overwhelming probability, the
[1]
random solution generated by Algorithm 3 is a k-exchange of Xt for some
k ≤ nγ ≤ n1/4 . This is again implied by the fact that with an overwhelming
γ
probability, at most nγ/2 low edges are chosen within the first T := n − 3n4
γ
steps in Algorithm 3, since the best case nγ/2 + 3n4 is still smaller than nγ .
By Claim 3, for any k ≤ nγ/2 and any m ≤ T, Algorithm 3 chooses
high edges with a probability at least 1 − 12/n at step m if at most k edges
have been chosen before step m, since there exist at least n − m − 3k ≥
3nγ
− 3nγ/2 ≥ 3 admissible high edges at step m.
4
19

Let P denote the probability of the random event that at most nγ/2 low
edges are chosen within T steps, and Q the probability of the random event
that at least nγ/2 + 1 low edges are chosen within the same T steps. Then
P = 1 − Q. We shall bound Q from above, which will give a lower bound for
P.
Let E be the random event that at least nγ/2 + 1 low edges are chosen
within T steps. Then Q = P[E]. For each l = 1, . . . , nγ/2 + 1, we define
a random variable vl denoting the first step m ≤ T such that l low edges
are chosen within m steps. Obviously, E implies the random event E1 that
v1 < v2 < · · · < vnγ/2 +1 ≤ T. Thus, Q ≤ P[E1 ], and P ≥ 1 − P[E1 ].
Observe that
X
P[E1 ] =
P[v1 = a1 , . . . , vnγ/2 +1 = anγ/2 +1 ],
a1 <a2 <···<anγ/2 +1 ≤T

and v1 = a1 , . . . , vnγ/2 +1 = anγ/2 +1 is equivalent to the random event that
before step a1 only high edges are chosen, that at any step between al and
al+1 only high edges are chosen for any l with 1 ≤ l ≤ nγ/2 , and that at steps
a1 , . . . , anγ/2 +1 only low edges are chosen. Thus, we have by Claim 3 that
P[v1 = a1 , . . . , vnγ/2 +1 = anγ/2 +1 ] ≤

12 nγ/2 +1
,
n

since at each step al , there exists at least one admissible high edge and we
do not care about what happens
 after step vnγ/2 +1 .
There are at most nγ/2T +1 different combinations for a1 < a2 < · · · <
 nγ/2 +1
.
anγ/2 +1 . Therefore, P ≥ 1 − P[E1 ] ≥ 1 − nγ/2T +1 12
n
γ/1
By Stirling’s formula, and observing that n + 1 ∈ o(T ), T ∈ Θ(n), we
 nγ/2 +1
γ/2
γ/2
have nγ/2T +1 12
= e−ω(n ) . Hence, P ≥ 1 − e−ω(n ) is overwhelmn
ingly large.
5. Main results
We shall now analyze the stochastic runtime of our two different random
solution generation methods for two classes of TSP instances that have been
well studied in the literature.

20

5.1. Stochastic runtime analysis for simple instances
We first consider a class of simple TSP instances that is defined by the
following distance function d : E → R on a graph with n vertices.


1 if {i, j} = {i, i + 1} for each i = 1, 2, . . . , n − 1,
d({i, j}) = 1 if {i, j} = {n, 1},
(10)


n otherwise.
Obviously, TSP instances with this distance function have a unique optimal
solution s∗ = (1, 2, . . . , n) (in the sense of the underlying Hamiltonian cycle),
and s∗ has a cost of n. The cost of an arbitrary feasible solution s equals
k + (n − k) · n, where k is the number of edges ∈ s that are also in s∗ . We
shall refer to these instances as G1 in the sequel.
The class G1 has been used in [19] and [42] for analyzing the expected
runtime of variants of MMAS. Zhou [19] proved that the (1 + 1) MMAA
algorithm has an expected runtime of O(n6 + n lnρ n ) on G1 in the case of nonvisibility (i.e., without the greedy distance information in the sampling), and
has an expected runtime of O(n5 + n lnρ n ) in the case of visibility (i.e., with
considering the greedy distance information in the sampling). Kötzing et al
[42] continued the study in [19]. They investigated the expected runtime of
(1 + 1) MMAA and its variant MMAS∗Arb on G1 and other TSP instances
on which both (1 + 1) MMAA and MMAS∗Arb have exponential expected
runtime. MMAS∗Arb differs with (1 + 1) MMAA only in the random solution
generation. MMAS∗Arb uses Algorithm 3 as its random solution generation
method, while (1 + 1) MMAA used Algorithm 2. Kötzing et al [42] proved
that MMAS∗Arb has an expected runtime of O(n3 ln n + n lnρ n ) on G1 .
Theorem 1 shows a stochastic runtime of O(n6+ ) for the CE variant with
the add-on, i.e., Algorithm 1 with max-min calibration (6), the vertex-based
random solution generation, and a stochastic runtime of O(n4+ ) for the
edge-based random solution generation. These results are comparable with
the above known expected runtimes. Although we are not able to get strictly
superior runtimes, our results are actually stronger and more informative.
Theorem 1 (Stochastic runtime of Algorithm 1 with max-min calibration
on G1 ). Assume that we set M = 1, ρ = 1, and use Algorithm 1 with the
1
, πmax = 1 − n1 . Then
max-min calibration (6) for the values πmin = n(n−2)
a) if we use the vertex-based random solution generation method (Algorithm 2), and take a sample size N ∈ Ω(n5+ ) for any constant  ∈
21

5

(0, 1), then with a probability at least 1 − e−Ω(N/n ) the optimal solution
s∗ can be found within n iterations;
b) if we use the edge-based random solution generation method (Algorithm
3), and take a sample size N ∈ Ω(n3+ ) for a constant  ∈ (0, 1), then
3
with a probability at least 1 − e−Ω(N/n ) , the optimal solution can be
found within n iterations.
Proof. We prove the Theorem by showing that the probability of the random
event that before the optimal solution is met, the number of edges shared by
the iteration-best and optimal solution strictly increases is overwhelmingly
large. This implies that the optimal solution is found within n iterations,
since the optimal solution has only n edges. Furthermore, the runtimes
presented in the Theorem hold. We only discuss the case of a), b) follows
with an almost identical argument.
[1]
By [19] (see also proof of Theorem 2), if Xt is not optimal, it has at
least either a 2-opt move or a 3-opt move. Note that for G1 , any k-opt
move of the iteration-best solution increases the number of its edges shared
with the optimal solution. By Claim 1, any 2-opt move is generated by
Algorithm 2 with probability Ω(n−3 ), and any 3-opt move is generated with
[1]
[1]
probability Ω(n−5 ). Thus, if Xt is not optimal, Xt+1 shares more edges with
[1]
the optimal solution than Xt with a probability at least 1 − (1 − n−5 )N =
5

1−e−Ω(N/n ) ∈ 1−e−Ω(n ) if N ∈ Ω(n5+ ) for any  > 0. Thus, this repeatedly
happens within polynomially many number of iterations with overwhelming
5
probability 1 − e−Ω(N/n ) . This completes the proof.
The stochastic runtimes of Theorem 1 are derived for a relatively large
sample size, namely N = Ω(n5+ ) and N = Ω(n3+ ). Actually, Theorem 1
may still hold for a smaller sample size. Theorem 2 partially asserts this.
It states that the total number of iterations required to reach the optimal
solution for both generation schemes may increase considerably if a smaller
sample size is used. However, the stochastic runtime does not increase. Interestingly, one can obtain a smaller stochastic runtime with a small sample
size for the edge-based random solution generation.
Theorem 2 (Stochastic runtime of Algorithm 1 on G1 for a small sample
size). Assume the conditions in Theorem 1, but set N ∈ Ω(n ) for any  ∈
(0, 1). Then:

22

a) For the vertex-based random solution generation, Algorithm 1 finds the
optimal solution s∗ within n6 iterations with a probability of 1 − e−Ω(N ) .
b) For the edge-based random solution generation, Algorithm 1 finds the
optimal solution s∗ within n3 ln n iterations with a probability of 1 −
e−Ω(N ) .
Proof of Theorem 2. The proof shares a similar idea with that of Theorem 1.
However, we consider here the random event that the number of edges shared
by the iteration-best and optimal solution does not decrease and strictly
increases enough times within a specified polynomial number of iterations.
For a), we shall consider the first n6 iterations. By Claim 2, the number
of edges shared by the iteration-best and optimal solution does not decrease
N
with a probability 1 − 1 − Ω(1)
= 1 − e−Ω(N ) (N ∈ Ω(n )). Therefore,
the number does not decrease within the first n6 iterations with probability
Qn6
−Ω(N )
) = 1−e−Ω(N ) . By Claim 1, for every consecutive n5 iterations,
t=0 (1−e
if the starting iteration-best solution is not optimal, then with probability
n5
1 − (1 − n−5 )N
= 1 − e−Ω(N ) , the number will strictly increase at least
once within these n5 iterations. Therefore, with overwhelming probability
1 − e−Ω(N ) , the optimal solution will be reached within the period of the first
n6 iterations, since there are n many consecutive n5 iterations within that
period.
b) can be proved by a similar way with a). We shall consider the first
3
n ln n iterations. By Claim 4, with probability 1 − (1 − Ω(1))N = 1 − e−Ω(N ) ,
the number of shared edges does not decrease in consecutive two iterations.
To complete the proof, we need an extra fact on 2, 3-exchanges.
Kötzing et al [42] showed for MMAS∗Arb that if the best solution s∗t found
so far has n−k edges from the optimal solution s∗ , then the probability of the
event that s∗t+1 has at least n−k+1 edges from s∗ , is in Ω(k/n3 ). We shall use a
different but simpler proof to show that this alsoTholds in our case of iteration[1]
best reinforcement. And with this fact, if |Xt
s∗ | = n−k for some 0 < k ≤
3
n, then with probability 1 − ((1 − k · n−3 )N )n /k = 1 − e−Ω(N ) , the number
of edges shared by the iteration-best solution and s∗ will strictly increase
at least once within the period [t, t + n3 /k]. This implies that s∗ is sampled
within the first n3 ln n iterations with overwhelming probability 1 − e−Ω(N ) ,
since n3 ln n iterations can be partitioned into n many consecutive phases
[0, n2 ), [n2 , n2 +n3 /(n−1)), [n2 +n3 /(n−1), n2 +n3 /(n−1)+n3 /(n−2)), . . . .
We now prove that fact.
23

[1]

We first show that when |Xt ∩ s∗ | = n − k with k > 0, then there exists a
[1]
2-opt move or a 3-opt move for Xt (see also [19] for a similar proof). Assume
[1]
that Xt contains exactly n − k edges from s∗ for some integer k > 0. Let
[1]
e∗ = {i, i + 1} be an edge in s∗ but not in Xt . Note that each node of the
[1]
graph is exactly incident to two edges of s∗ and Xt , respectively. Therefore
[1]
[1]
there exists an edge e0 ∈ Xt incident to i, an edge e00 ∈ Xt incident to i+1,
and e0 , e00 are not in s∗ . Figure 2 shows an example, where e0 is either {i, u} or
{i, v}, and e00 is either {i+1, w} or {i+1, y}. If e0 = {i, u} and e00 = {i+1, w} or

u

i

e

v

i+1

∗

w

y
e1

[1]

Figure 2: Demonstration of adding a new edge. The solid edges represent the cycle Xt .
[1]

if e0 = {i, v} and e00 = {i+1, y}, then there exists a 2-opt move of Xt which
removes e0 , e00 of distance n and adds e∗ and another edge (either {u, w} or
{v, y}) of distance at most n + 1 together. If e0 = {i, u}, e00 = {i+1, y} or
[1]
e0 = {i, v}, e00 = {i+1, w}, there is a 3-opt move of Xt which removes e0 , e00 ,
and an edge e1 ∈
/ s∗ , and adds edge e∗ and another two edges, this replacing
3 edges of distance n by 3 edges of distance at most 2n+1 together. Here,
[1]
[1]
observe the fact that adding e∗ to Xt and removing e0 , e00 from Xt results
[1]
in graph containing a cycle, and there must be an edge e1 ∈ Xt on that cycle
that does not belong to s∗ . We choose this edge as the edge e1 . Therefore,
[1]
for each e∗ of the k remaining edges in s∗ that are not in Xt , there exists a
[1]
2-opt or 3-opt move of Xt that adds e∗ .
By Claim 4, for any l ∈ O(1), the probability of producing an l-exchange
[1]
of the iteration-best solution Xt by Algorithm 3 in iteration t + 1 is Ω(1).
Since any two l-exchanges are produced with the same probability, the probability of producing a particular l-exchange in iteration t + 1 is Ω(1/nl ). As
[1]
a result, Algorithm 3 produces for each edge e∗ ∈ s∗ − Xt a 2-opt or 3-opt
[1]
move of Xt that adds edge e∗ with probability at least Ω(1/n3 ).
24

Note that the generation of a 2-exchange (or a 3-exchange) with two newly
added edges e2 , e3 by Algorithm 3 includes two mutually exclusive cases (3!
cases for a 3-exchange): e2 is chosen before e3 , or e3 is chosen before e2 . It
is not difficult to see that these two cases (3! cases for 3-exchange) have the
same probability. Therefore, the probability of the event that Algorithm 3
[1]
generates a 2-opt or 3-opt move of Xt that e∗ as one of the newly added edges
and selects e∗ before the other newly added edges, is bounded from below
[1]
by Ω(1/3!n3 ) = Ω(1/n3 ). Since Xt has k such e∗ and the corresponding k
[1]
events are also mutually exclusive, we obtain that the probability that Xt+1
[1]
[1]
has more edges from s∗ than Xt if Xt has exactly n − k edges from s∗ for
a constant k > 0 is Ω(k/n3 )
Corollary 1 further improves the stochastic runtime for an even smaller
sample size. It can be proved by an argument similar to the proof of Theorem
l
1, where we observe that (1 − (1 − p(n))ω(ln n) )n = 1 − n−ω(1) for any constant
l > 0 and probability p(n) ∈ Ω(1), and that 1 − e−ω(ln n) = 1 − n−ω(1) .
Corollary 1. Assume the conditions in Theorem 1, but let N ∈ ω(ln n).
Then:
a) For the vertex-based random solution generation, Algorithm 1 finds the
optimal solution s∗ within n6 iterations with a probability of 1 − n−ω(1) .
Particularly, if N = (ln n)2 , the runtime is n6 (ln n)2 with probability
1 − n−ω(1) .
b) For the edge-based random solution generation, Algorithm 1 finds the
optimal solution s∗ within n3 ln n iterations with a probability of 1 −
n−ω(1) . Particularly, if N = (ln n)2 , the runtime is n3 (ln n)3 with probability 1 − n−ω(1) .
Theorem 2 tells that, for any  ∈ (0, 1), a sample size of N ∈ Θ(n ) is
already sufficient for iteration-best reinforcement to efficiently find an optimal
solution of simple TSP instances with an overwhelming probability. Corollary
1 further shows that N ∈ ω(ln n) even leads to a better runtime with a slightly
smaller but still overwhelming probability. Theorem 3 below shows that with
an overwhelming probability, the runtime of iteration-best reinforcement will
be exponential if N ∈ O(ln n), even if the instances are as simple as those in
G1 .
25

1
ln n.
Theorem 3. Assume the conditions of Theorem 1, but set N < 220
1/200
−Ω(n
)
Then, with probability 1 − e
, Algorithm 1 with edge-based solution
1/300 )
generation does not find the optimal solution s∗ within eΘ(n
iterations.

Proof. We prove the Theorem by inspecting the probability of the random
event that, before the optimal solution is found, the cost of the iteration[1]
best solution Xt will oscillate for exponentially many iterations with an
overwhelming probability. We shall consider this in the last stages of the
optimization process.
Let T0 be the first iteration which samples a solution containing at least
n − n1/4 + n1/5 edges from the optimal solution. We show that with an
overwhelming probability, the number of common edges in the iteration-best
and optimal solution will drop below n − n1/4 + n1/5 and the optimal solution
is not sampled before that. This will imply the conclusion of Theorem 3,
since, with an overwhelming probability, this phenomenon can repeatedly
occur exponentially many times before optimal solution is found.
To that end, we need to show the following:
[1]

1) For any 1/4 > δ > 0, if Xt contains at least n − nδ edges from the
optimal solution, then with a probability O( √1n ), the random solution
generated by Algorithm 3 will contain more edges from the optimal
[1]
solution than Xt in iteration t + 1;
[1]

2) For any 1/4 > δ > 0, if Xt contains at least n − nδ edges from
the optimal solution, then with a probability Ω(1) (at least e−5 ), the
random solution generated by Algorithm 3 will contain fewer edges
[1]
from the optimal solution than Xt in iteration t + 1.
However, we first use these two facts and show them afterwards.
1/20
[1]
By Claim 6, with probability 1−e−ω(n ) , XT0 contains at most n−n1/4 +
n1/5 + n1/10 edges from the optimal solution, since the random event that
the number of common edges from the iteration-best and optimal solution
increases more than n1/10 in one iteration implies an occurrence of a Ω(n1/10 )1/200 )
exchange. Similarly, by Claim 6 again, with probability 1 − e−ω(n
, the
iteration-best solution contains k ∈ [n − n1/4 + n1/5 − n1/6+1/100 , n − n1/4 +
n1/5 + n1/10 + n1/6+1/100 ] edges from the optimal solution in each iteration
t ∈ [T0 , T0 + n1/6 ]. This means that the optimal solution is not found in the
period [T0 , T0 + n1/6 ] with an overwhelming probability. With the help of
1) and 2), we are now to show that within this period, the number of edges
26

shared by the iteration-best and optimal solution is significantly reduced with
an overwhelming probability. This will complete the proof.
To facilitate our discussion, we call an iteration a successful iteration
if its iteration-best solution contains more edges from the optimal solution
than the last iteration-best solution, and an iteration a failure iteration if its
iteration-best solution contains fewer edges from the optimal solution than
the last iteration-best solution.
By 1) and the subsequent discussion, the expected number of success1
n
), since N < 220
ln n. Thus,
ful iterations within [T0 , T0 + n1/6 ] is O( nln1/3
1/6
−Ω(n )
by the Chernoff bound, with probability 1 − e
, at most n1/100 successful iterations can occur within [T0 , T0 + n1/6 ]. By 2) and the subsequent discussion, the expected number of failure iterations in [T0 , T0 + n1/6 ]
1
1
1
ln n. By the Chernoff bound, it happens that
is Ω(n 6 − 44 ), since N < 220
−Ω(n1/6 )
with probability 1 − e
, at least n1/7 failure iterations will occur in
[T0 , T0 + n1/6 ]. Since a successful iteration can add at most n1/100 edges from
1/200 )
the optimal solution with probability 1 − e−ω(n
, it totally adds at most
1/100
1/100
1/50
n
×n
=n
edges from the optimal solution to the iteration-best
1/200 )
solution within [T0 , T0 +n1/6 ] with probability 1−e−ω(n
. Note that within
1/6
−Ω(n1/6 )
1/7
[T0 , T0 + n ], with probability 1 − e
, at least n × 1 = n1/7 “good”
edges are removed from the iteration-best solution. Therefore, with over1/200 )
[1]
whelming probability 1 − e−Ω(n
, XT0 +n1/6 will contain at most
n − n1/4 + n1/5 + n1/10 − n1/7 + n1/50 < n − n1/4 + n1/5
[1]

edges from the optimal solution, since XT0 contains at most n − n1/4 + n1/5 +
1/20
n1/10 iterations with probability 1 − e−Ω(n ) . As a result, with probability
1/200 )
1−e−Ω(n
, the number of common edges in the iteration-best and optimal
solution will again be smaller than n − n1/4 + n1/5 in some iteration after T0 ,
and the optimal solution is not found before that. And this will repeatedly
1/300 )
1/200 )
happen eΘ(n
times with probability 1 − e−Ω(n
.
To finish the proof, we now formallyprove 1) and 2). We first consider 2).
n
By taking
√ k = 2 and considering the 2 2-exchanges that happen in the first
n − 3 n steps in the proof of Claim 4, one can show a tighter probability
[1]
lower bound e15 for producing 2-exchanges of Xt by Algorithm 3. Here, √
we
observe that the probability of choosing a high edge at a step before n − 3 n
is at least 1 − 3/(n − 2), see the proof of Claim 3.
Note that if 2-exchanges deleting 2 edges from the optimal solution happen N times in an iteration, then the iteration will be a failure iteration.
27

By the above and the fact that any two k-exchanges happen with the same
probability, a failure iteration then occurs with a probability at least
1
e5

n−nδ
2
n
2

 !N
≥

1
e5

n−nδ
2
n
2

1
 ! 220
ln n

∈ Ω(n−1/44 ),

1
where δ ∈ (0, 1/4) and N < 220
ln n. This asserts 2).
[1]
1) follows with a similar discussion. Since Xt is assumed to contain
at least n − nδ edges from the optimal solution for some δ ∈ (0, 1/4), and
since Ω(nδ )-exchanges happen with an overwhelmingly small probability, we
need to consider only O(nδ )-exchanges when we estimate the probability of a
successful iteration. For each k ∈ Ω(nδ ), the proportion of failure k-exchanges
is bounded from below by

n−nδ
2knδ
−1/2
k
= e− n + o(1) ≥ e−2n
+ o(1),
n
k

since 0 < δ < 1/4, and k-exchanges removing k edges shared by the iterationbest and optimal solution are not “successful” k-exchanges. Since for any
k ∈ Ω(nδ ), any two k-exchanges happen with the same probability, and
since the sum of the probabilities of successful and failure k-exchanges is
smaller than 1, we conclude that successful O(nδ )-exchanges happen with a
−1/2
probability smaller than 1−e2n
∈ O( √1n ). Therefore, a successful iteration
√ n ) since N < 1 ln n.
happens with a probability 1−(1−O( √1n ))N ∈ O( ln
220
n
Theorem 3 generalizes the finding of [20] to simple TSP instances. It
formally states that for ρ = 1, N ∈ Ω(ln n) is necessary to efficiently find
an optimal solution to TSP. By Theorem 3, Theorem 1, Theorem 2 and its
Corollary 1, we have clearly analyzed the impact of the size of N on the
resulting stochastic runtime for the simple TSP instances in the case of that
ρ = 1. N ∈ ω(ln n) is sufficient to find the optimal solution in a stochastically
polynomial runtime, and the degree of the polynomial may increase with N ,
but the probability guaranteeing the runtime is also increasing with N .
5.2. Stochastic runtime analysis for grid instances
Now, we consider more general TSP instances. Herein, the n vertices are
positioned on an m × m grid for some integer m ∈ N+ . The vertices are
positioned in a way that no three of them are collinear. Figure 3 gives an
28

5
4
3
2
1

0

1

2

3

4

5

Figure 3: A grid instance

example of such an instance where m = 5 and n = 8. The weight of an edge
{l, k} ∈ E in this case is defined as the usual Euclidean distance d(l, k)
between vertex l and vertex k for every l, k = 1, . . . , n. In this section, we
shall refer to these TSP instances as grid instances.
Grid instances have been studied in [43] and [30]. Sutton and Neumann
[43] investigated the expected runtime of (1+1) EA and RLS for these instances. As a continuation of [43], Sutton et al [30] further proved that
the more extensive algorithm (µ + λ) EA finds an optimal solution for the
instances expectedly in
O((µ/λ)n3 m5 +nm5 +(µ/λ)n4k (2k−1)!)
iterations if every of the λ selected parents is mutated by taking a random
number of consecutive 2-exchange moves, and expectedly in
O((µ/λ)n3 m5 +nm5 +(µ/λ)n2k (k−1)!)
iterations with a mixed mutation operator, where k denotes the number of
vertices that are not on the boundary of the convex hull of V. Sutton et
al [30] also studied general Euclidean TSP instances (without collinearity)
and showed similar results in terms of the maximum distance value dmax ,
the minimum distance value dmin , k and the minimum angle in the triangles
formed by the vertices.
Before we present our stochastic runtime, we summarize some structural
properties of grid instances (some just follow from properties of general Euclidean instances). We say that two different edges {i, j} and {k, l} intersect
with each other if there exists a point p such that p ∈
/ {i, j, k, l} locates
on both of the two edges, see, e.g., Figure 4a. We say that a solution is
29

intersection-free if the corresponding Hamiltonian cycle does not contain intersections, see, e.g., Figure 4b.

j
k
p

l

j
l

i
(a) intersection

k

i

(b) intersection free

Figure 4: Example for intersections

Obviously, the triangle inequality [49] holds for grid instances. Therefore,
removing an intersection by a (unique) 2-exchange move in a solution strictly
reduces the total traveling cost, see Figure 4a. Lemma 1 states the well known
fact that an optimal solution of grid instances is intersection-free.
Lemma 1. Optimal solutions of grid instances are intersection-free.

i

l
p

k

j

Figure 5: Example for a 2-opt move

We now restrict 2-opt moves to 2-exchange moves that remove an intersection. For example, removing edges {i, j}, {k, l} in Figure 5 and adding
new edges {i, l}, {k, j} form such a 2-opt move. Lemma 2 below says that
for grid instances, removing one intersection may reduce the total traveling
cost Ω(m−4 ) if it is applicable. We omit the simple proof here. Interested
readers may refer to [30] for a proof.
Lemma 2. If a feasible solution to a grid instance contains intersections,
then removing the intersection can reduce the total traveling cost Ω(m−4 ).
30

The convex hull Y(V ) of the vertex set V is the smallest convex set in R2
that contains V . Its boundary is a convex polygon spanned by some vertices
with possibly other vertices in the interior of that polygon. Let V b denote
the set of vertices on the boundary of Y(V ). Figure 6 illustrates this.

Figure 6: Example of a convex hull

Quintas and Supnick [50] proved that if a solution s is intersection-free,
then the solution respects the hull-order, i.e., any two vertices in the subsequence of s induced by the boundary (the outer polygon) of Y(V ) are
consecutive in s if and only if they are consecutive on the boundary of Y(V ).
Therefore, if V b = V, i.e., all of the vertices are on the convex hull, then
every intersection-free solution is optimal.
Theorem 4 below analyzes the stochastic runtime of Algorithm 1 for grid
instances for the case that V = V b . It states that the stochastic runtime is
O(n4 ·m5+ ) for the vertex-based random solution generation, and O(n3 ·m5+ )
for the edge-based random solution generation. Corollary 2 further improves
the runtime by sacrificing the probability guarantee. These stochastic runtimes are close to the expected runtime O(n3 ·m5 ) for RLS reported by Sutton
et al [43] and [30].
Theorem 4. Consider a TSP instance with n vertices located on an m × m
grid such that no three of them are collinear. Assume that V b = V , i.e., every
vertex in V is on the convex hull V b , that we apply the max-min calibration
1
, ρ = 1, M = 1 and N ∈ Ω(m ) for some
(6) with πmax = 1 − n1 , πmin = n(n−2)
constant  > 0. Then:
a) With an overwhelming probability of 1 − e−Ω(N ) , Algorithm 1 finds the
optimal solution within at most n4 · m5 iterations with the vertex-based
random solution generation.
b) With an overwhelming probability of 1 − e−Ω(N ) , Algorithm 1 finds an
optimal solution within at most n3 · m5 iterations with edge-based random solution generation.
31

Proof of Theorem 4. Note that under the conditions of Theorem 4, every
intersection free solution is optimal. By Lemma 2, we know that a 2-opt move
reduces the total traveling cost by Ω(m−4 ). Therefore, n·m5 consecutive 2-opt
moves turn a feasible solution into an optimal one, since the worst solution in
this case has a total traveling cost smaller than n·m and the optimal solution
has total traveling cost larger than n. Notice also that m ≥ n/2, since the n
vertices are positioned on the m × m grid and no three of them are collinear.
With these facts, we prove the Theorem by a similar argument to the one
used in the proof of Theorem 2.
Again, we consider the random event that the cost of the iteration best
solution does not increase within a specified period of polynomially many
iterations and strictly decreases sufficiently many times within that period.
For a), we consider the first n4 m5 iterations. For b), we consider the first
n3 m5 iterations.
4 5
For a) : By Claim 2, with probability (1 − (1 − Ω(1))N )n m = 1 − e−Ω(N ) ,
the cost of the iteration-best solution does not increase within n4 m5 iterations. By Claim 1, for a phase consisting of consecutive n3 iterations, with
3
probability 1 − (1 − n−3 )N ·n = 1 − e−Ω(N ) , in at least one iteration of that
phase an intersection is removed from the iteration-best solution, provided
the phase starts with an iteration-best solution containing at least one intersection. Since the first n4 m5 iterations can have nm5 such phases, a)
follows.
b) follows with an almost identical discussion. We therefore omit the
proof.
Corollary 2. Consider a TSP instance with n vertices located on an m × m
grid such that no three of them are collinear. Assume that V b = V , i.e., every
vertex in V is on the convex hull V b , that we apply the max-min calibration
1
, ρ = 1, M = 1 and N ∈ ω(ln m). Then:
(6) with πmax = 1 − n1 , πmin = n(n−2)
a) With probability 1 − m−ω(1) , Algorithm 1 finds the optimal solution
within at most n4 · m5 iterations with the vertex-based random solution generation.
b) With probability 1−m−ω(1) , Algorithm 1 finds an optimal solution within
at most n3 · m5 iterations with the edge-based random solution generation.
Now, we consider the more interesting case that |V | − |V b | = k ∈ O(1),
i.e., k vertices are not on the convex hull. Note that we can turn an arbitrary
32

intersection-free solution to an optimal solution only by rearranging the positions of those k interior points in that solution, and this requires at most k
consecutive jump moves (see [30] for a proof). A jump move δi,j transforms a
solution into another solution by shifting positions i, j as follows. Solution s
is transformed into solution δi,j (s) by moving the vertex at position i into position j while vertices at positions between i and j are shifted appropriately,
e.g.,
δ2,5 (i1 , i2 , i3 , i4 , i5 , i6 , i7 ) = (i1 , i3 , i4 , i5 , i2 , i6 , i7 ) and
δ5,2 (i1 , i2 , i3 , i4 , i5 , i6 , i7 ) = (i1 , i5 , i2 , i3 , i4 , i6 , i7 ).
It is not difficult to see that a jump move δi,j can be simulated by either a
2-exchange move (in the case that |i−j| = 1) or a 3-exchange move (in all
other cases). Therefore, we can actually turn an intersection-free solution
into an optimal one by a sequence of at most k consecutive 2-exchange or
3-exchange moves. Furthermore, a sequence of k consecutive 2-exchange or
3-exchange moves can be simulated by a κ-exchange move with an integer
κ ≤ 3k. This means that any intersection-free solution can be turned into an
optimal solution by a κ-exchange move with κ ≤ 3k. We shall call such a κexchange move in the sequel a 3k-opt move, although κ may be smaller than
1
) by
3k. Recall that a 3k-opt move is produced with a probability of Ω( n6k−1
1
Algorithm 2 (see Claim 1), and with a probability of Ω( n3k ) by Algorithm
3 (see Lemma 6 of [42], or Claim 4) in any of the N independent draws in
[1]
iteration t, if Xt−1 is intersection-free and not optimal. As a result, we obtain
by a similar proof as above Theorem 5 below.
Theorem 5. Consider a TSP instance with n vertices located on an m × m
grid such that no three of them are collinear. Assume that |V | − |V b | = k ∈
O(1) (k vertices are not on the convex hull V b ), that we apply the max-min
1
calibration 6 with πmax = 1 − n1 , πmin = n(n−2)
, and set ρ = 1, M = 1, for
some constant  > 0. Then:
a) If we set N ∈ Ω(n3 · m ), then with an overwhelming probability of
3
1 − e−Ω(N/n ) , Algorithm 1 finds an optimal solution within at most n ·
m5 + n6k−4 iterations with the vertex-based random solution generation;
b) If we set N ∈ Ω(n2 · m ), then with an overwhelming probability of
2
1 − e−Ω(N/n ) , Algorithm 1 finds an optimal solution within at most
n·m5 +n3k−2 iterations with the edge-based random solution generation.

33

Proof of Theorem 5. We only prove a). b) can be derived by a very similar
argument. We define two random events as following:
[1]

[1]

E1 : for each t ≤ n · m5 + n6k−4 , f (Xt−1 ) ≥ f (Xt );
[1]

E2 : for each t ≤ n · m5 + n6k−4 , if Xt−1 is not intersection-free, then a 2-opt
move happens in iteration t.
By a similar argument as the one for Theorem 4, we obtain that P[E1 ∩ E2 ] ≥
3
1−e−Ω(N/n ) . Let η be a random variable denoting the number of iterations for
[1]
which Xt−1 is intersection-free. Notice that, conditioned on E1 ∩E2 , η ≥ n·m5
implies that an optimal solution occurs within n · m5 + n6k−4 iterations.
Conditioned on E1 ∩ E2 and η < n · m5 , there are at least Ω(n6k−4 ) itera[1]
tions in which Xt−1 is intersection-free, since each Xt−1 is either intersection[1]
free or not intersection-free. Note also that in each iteration in which Xt−1
[1]
intersection-free and not optimal, a 3k-opt move that turns Xt−1 into an op1
))N . This
timal solution happens with probability of at least 1 − (1 − Ω( n6k−1
[1]
means for any fixed t ∈ N, if Xt−1 is intersection-free, then the probability of
[1]
1
the event that Xt is optimal is bounded from below by 1 − (1 − Ω( n6k−1
))N .
Therefore, for any fixed Ω(n6k−4 ) iterations in which the iteration-best solu[1]
tion Xt−1 is intersection-free and not optimal, the probability of the event
[1]
that the corresponding Ω(n6k−4 ) Xt ’s are still not optimal, is bounded from
6k−4
3
1
above by (1 − Ω( n6k−1
))N ·n
= e−Ω(N/n ) . This means that, conditioned on
E1 ∩ E2 and η < n · m5 , an optimal solution occurs within n · m5 + n6k−4
3
iterations with a probability of 1 − e−Ω(N/n ) .
As a result, an optimal solution occurs within the first n · m5 + n6k−4
3
iterations with a probability of 1 − e−Ω(N/n ) .
Theorem 5 shows a stochastic runtime of n3 m5+ +n6k−1 m for Algorithm
1 equipped with the vertex-based solution generation, and a stochastic runtime of n3 m5+ + n3k m for Algorithm 1 equipped with edge-based solution
generation, in the case of that |V | − |V b | = k ∈ O(1). This is much better
than the expected runtime
O(µ · n3 m5 +nm5 +µ · n4k (2k−1)!)
for (µ+λ) EA with sequential 2-opt mutations reported by Sutton et al [30].
However, we are not able to analyze the stochastic runtime in the case that
34

k ∈ ω(1), since k ∈ ω(1) interior points may require super-polynomially
many iterations to turn an intersection-free solution into an optimal solution
when a polynomial sample size is used.
6. Conclusion
We have analyzed the stochastic runtime of a CE algorithm on two classes
of TSP instances under two different random solution generation methods.
The stochastic runtimes are comparable with corresponding expected runtimes reported in the literature.
Our results show that the edge-based random solution generation method
makes the algorithm more efficient for TSP instances in most cases. Moreover, N ∈ Ω(ln n) is necessary for efficiently finding an optimal solution with
iteration-best reinforcement. For simple instances, N ∈ ω(ln n) is sufficient
to efficiently find an optimal solution with an overwhelming probability, and
N ∈ O(ln n) results in an exponential runtime with an overwhelming probability. However, for more difficult instances, one may need to use a relatively
large sample size.
Our stochastic runtimes are better than the expected runtimes of the
(µ + λ) EA on the grid instances. The EA randomly changes local structures
of some of its current solutions by a Poisson distributed number of consecutive 2-exchange moves in every iteration, while our algorithm refrains from
local operations on current solutions and only refreshes solutions by sampling from an evolving distribution. The solution reproducing mechanism in
the EA stays the same throughout the optimization, only the current solutions in every iteration vary. However, the solution reproducing mechanism
(sampling distribution) of our algorithm also evolves. This is the essential
difference of MBS with traditional EAs. The comparison of our results with
the expected runtimes in [30] therefore show that using a self-adaptive dynamic solution reproducing mechanism is helpful (in efficiently finding an
optimal solution) when the search space becomes rugged. The stochastic
runtimes in Theorem 4 are only valid for instances with a bounded number
of interior points. In the future, it should be interesting to analyze the case
that |V | − |V b | ∈ ω(1). This might also give more insight to the problem of
RP v.s. P [51].
Our analysis is actually a kind of worst-case analysis, which is rather
pessimistic. We analyze the optimization progress by only checking some very
particular random events. This may not only underestimate the probability
35

of finding an optimal solution with our algorithm, but also overestimate the
required number of iterations. In the future, it should be of great interest to
consider a smoothed runtime analysis over an -neighborhood of the n nodes
in the real plane as has been done for the Simplex method by Spielman and
Teng in their famous paper [52].
Acknowledgment
We thank the anonymous reviewers for their numerous useful suggestions
on improving the scientific quality and English presentation of this article.
References
[1] R. Y. Rubinstein, D. P. Kroese, The cross-entropy method: a unified
approach to combinatorial optimization, Monte-Carlo simulation and
machine learning, Springer Science & Business Media, 2004.
[2] R. Y. Rubinstein, Optimization of computer simulation models with rare
events, European Journal of Operational Research 99 (1) (1997) 89–112.
[3] R. Y. Rubinstein, The cross-entropy method for combinatorial and continuous optimization, Methodology and computing in applied probability 1 (2) (1999) 127–190.
[4] M. Dorigo, T. Stützle, Ant colony optimization, Cambridge, Massachusetts: A Bradford Book, MIT Press, 2004.
[5] M. Hauschild, M. Pelikan, An introduction and survey of estimation
of distribution algorithms, Swarm and Evolutionary Computation 1 (3)
(2011) 111–128.
[6] M. Zlochin, M. Birattari, N. Meuleau, M. Dorigo, Model-based search
for combinatorial optimization: A critical survey, Annals of Operations
Research 131 (1-4) (2004) 373–395.
[7] D. Whitley, A genetic algorithm tutorial, Statistics and computing 4 (2)
(1994) 65–85.
[8] H. R. Lourenço, O. C. Martin, T. Stützle, Iterated local search, Springer,
2003.
36

[9] Z. Wu, Model-based heuristics for combinatorial optimization: a mathematical study of their asymptotic behavior, Ph.D. thesis, Institut für
Angewandte Stochastik und Operations Research (IASOR), Technical
University of Clausthal (2015).
[10] T. Stützle, H. H. Hoos, MAX-MIN ant system, Journal of Future Generation Computer Systems 16 (2000) 889–914.
[11] S. Droste, T. Jansen, I. Wegener, On the analysis of the (1+1) evolutionary algorithm, Theoretical Computer Science 276 (1-2) (2002) 51–81.
[12] J. He, X. Yao, Drift analysis and average time complexity of evolutionary
algorithms, Artificial Intelligence 127 (1) (2001) 57–85.
[13] F. Neumann, C. Witt, Runtime analysis of a simple ant colony optimization algorithm, Tech. rep., Departmant of Computer Science, University
of Dortmund, Germany (2006).
[14] C. Witt, Runtime analysis of the (µ +1) ea on simple pseudo-boolean
functions, Evolutionary Computation 14 (1) (2006) 65–86.
[15] F. Neumann, C. Witt, Runtime analysis of a simple ant colony optimization algorithm, Algorithmica 54 (2) (2009) 243–255.
[16] B. Doerr, F. Neumann, D. Sudholt, C. Witt, Runtime analysis of the 1ant ant colony optimizer, Theoretical Computer Science 412 (17) (2011)
1629–1644.
[17] W. J. Gutjahr, G. Sebastiani, Runtime analysis of ant colony optimization with best-so-far reinforcement, Methodology & Computing in Applied Probability 10 (3) (2008) 409–433.
[18] Y. Zhou, J. He, A runtime analysis of evolutionary algorithms for constrained optimization problems, IEEE Transactions on Evolutionary
Computation 11 (5) (2007) 608–619.
[19] Y. Zhou, Runtime analysis of an ant colony optimization algorithm for
tsp instances, Evolutionary Computation IEEE Transactions on 13 (5)
(2009) 1083–1092.

37

[20] F. Neumann, D. Sudholt, C. Witt, A few ants are enough:aco with
iteration-best update, in: Genetic and Evolutionary Computation Conference, GECCO 2010, Proceedings, Portland, Oregon, Usa, July, 2010,
pp. 63–70.
[21] P. S. Oliverto, C. Witt, Improved time complexity analysis of the simple
genetic algorithm, Theoretical Computer Science 605 (15) (2015) 21–41.
[22] D. Sudholt, C. Thyssen, Runtime analysis of ant colony optimization for
shortest path problems, Journal of Discrete Algorithms 10 (10) (2012)
165–180.
[23] A. Lissovoi, C. Witt, Runtime analysis of ant colony optimization on dynamic shortest path problems, Theoretical Computer Science 561 (2015)
73–85.
[24] Y. Chen, X. Zou, Runtime analysis of a multi-objective evolutionary algorithm for obtaining finite approximations of pareto fronts, Information
Sciences 262 (2014) 62–77.
[25] D. Sudholt, C. Witt, Update strength in edas and aco: How to avoid
genetic drift, in: Genetic and Evolutionary Computation Conference,
2016, pp. 61–68.
[26] D. H. Wolpert, W. G. Macready, No free lunch theorems for optimization, IEEE Transactions on Evolutionary Computation 1 (1) (1997) 67–
82.
[27] F. Neumann, D. Sudholt, C. Witt, Analysis of different mmas aco algorithms on unimodal functions and plateaus, Swarm Intelligence 3 (2009)
35–68.
[28] F. Neumann, I. Wegener, Randomized local search, evolutionary algorithms, and the minimum spanning tree problem, Theoretical Computer
Science 378 (2007) 32–40.
[29] J. Reichel, M. Skutella, Evolutionary algorithms and matroid optimization problems, Algorithmica 57 (1) (2010) 187–206.
[30] A. M. Sutton, F. Neumann, S. Nallaperuma, Parameterized runtime
analyses of evolutionary algorithms for the planar euclidean traveling
salesperson problem, Evolutionary Computation 22 (4) (2014) 595–628.
38

[31] A. M. Sutton, J. Day, F. Neumann, A parameterized runtime analysis
of evolutionary algorithms for max-2-sat, in: Conference on Genetic &
Evolutionary Computation, 2012, pp. 433–440.
[32] Y. Zhou, X. Lai, K. Li, Approximation and parameterized runtime analysis of evolutionary algorithms for the maximum cut problem., IEEE
Transactions on Cybernetics 45 (8) (2015) 1491–1498.
[33] Z. Wu, M. Kolonko, Asymptotic properties of a generalized cross entropy
optimization algorithm, IEEE Transactions on Evolutionary Computation 18 (5) (2014) 658 – 673.
[34] Z. Wu, M. Kolonko, R. H. Möhring, Stochastic runtime analysis of the
cross entropy algorithm, IEEE Transactions on Evolutionary Computation, DOI: 10.1109/TEVC.2017.2667713.
[35] M. Held, R. M. Karp, A dynamic programming approach to sequencing
problems, Journal of the Society for Industrial and Applied Mathematics
10 (1) (1962) 196–210.
[36] N. Christofides, Worst-case analysis of a new heuristic for the travelling
salesman problem, Tech. rep., Graduate School of Industrial Administration, CMU (1976).
[37] M. T. Goodrich, R. Tamassia, Algorithm Design and Applications, Wiley, 2015.
[38] S. Arora, Polynomial time approximation schemes for Euclidean traveling salesman and other geometric problems, Journal of the ACM 45 (5)
(1988) 753–782.
[39] J. S. B. Mitchell, A constant-factor approximation algorithm for tsp with
pairwise-disjoint connected neighborhoods in the plane, in: TwentySixth Symposium on Computational Geometry, 2010, pp. 183–191.
[40] S. Lin, B. W. Kernighan, An effective heuristic algorithm for the
traveling-salesman problem, Operations Research 21 (2) (1973) 498–516.
[41] P. T. D. Boer, D. P. Kroese, S. Mannor, R. Y. Rubinstein, A tutorial
on the cross-entropy method, Annals of Operations Research 134 (1)
(2005) 19–67.
39

[42] T. Kötzing, F. Neumann, H. Röglin, C. Witt, Theoretical analysis of two
aco approaches for the traveling salesman problem, Swarm Intelligence
6 (1) (2012) 1–21.
[43] A. M. Sutton, F. Neumann, A parameterized runtime analysis of evolutionary algorithms for the euclidean traveling salesperson problem, in:
Proceedings of the Twenty-Sixth Conference on Artificial Intelligence
(AAAI’12), AAAI press, 2012, pp. 1105–1111.
[44] A. Costa, O. D. Jones, D. Kroese, Convergence properties of the crossentropy method for discrete optimization, Operations Research Letters
35 (5) (2007) 573–580.
[45] Z. Wu, M. Kolonko, Absorption in model-based search algorithms for
combinatorial optimization, in: Evolutionary Computation (CEC), 2014
IEEE Congress on, IEEE, 2014, pp. 1744–1751.
[46] M. Thomas, Machine learning, New Delhi: McGraw Hill Education India, 1997.
[47] H. Asoh, H. Mühlenbein, On the mean convergence time of evolutionary
algorithms without selection and mutation, in: Parallel Problem Solving
from Nature—PPSN III, Springer, 1994, pp. 88–97.
[48] M. Pirlot, General local search methods, European Journal of Operational Research 92 (3) (1996) 493–511.
[49] M. A. Khamsi, W. A. Kirk, An introduction to metric spaces and fixed
point theory, John Wiley,, 2001.
[50] L. V. Quintas, F. Supnick, On some properties of shortest Hamiltonian
circuits, American Mathematical Monthly 72 (9) (1965) 977–980.
[51] W. Gasarch, Classifying problems into complexity classes, Advances in
Computers 95 (2015) 239–292.
[52] D. A. Spielman, S. H. Teng, Smoothed analysis of algorithms: Why the
simplex algorithm usually takes polynomial time, Journal of the Acm
51 (3) (2004) 385–463.

40

