## The Traveling Salesman Problem in Theory & Practice

### Lecture 8: Lin-Kernighan and Beyond 11 March 2014

David S. Johnson [email protected]

http://davidsjohnson.net

Seeley Mudd 523, Tuesdays and Fridays

## 4-Opt

Seems to be relatively straightforward to generalize from 3-opt.

– – – Just need to now consider possibilities for t 7 and t 8 .

Can further exploit the Partial Sum Theorem.

Can consider all possibilities for t5, but need to make sure that t7 breaks up the short cycle (or create a new one).

An alternative View: Maintain the path

t 1 t 3 t 4 t 2 t 4

An alternative View: Maintain the path

t 1 t 3 t 2 t 4

An alternative View: Maintain the path, at least eventually…

t 1 t 4 t 3 t 5 t 6 t 2

t

An alternative View: Maintain the path, at least eventually…

1 t 4 t t 6 2

## Double Bridge Move

Cannot be produced by any sequential move.

Best one can be found in O(N 2 ) time.

2

### )

[Glover, “Finding a best traveling salesman 4-opt move in the same time as a best 2-opt move,” J. Heuristics 2 (1996), 169-179] Simpler algorithm due to Johnson & McGeoch, 2002

### Normal Form for Double Bridge Move

j j+1 p p+1

Smallest index of the first endpoint of a deleted edge

qp+1 q i+1 i

2 1

N N-1

C[h,k] = Cost of move that deletes {h,h+1} and {k,k+1} and adds {h,k+1} and {h+1,k}

## The Functions to be Computed

• • f(p,j) = min {C[p,q] : j < q ≤ N}, 2 ≤ p < j ≤ N .

g(i,j) = C[i,j] + min {f[p,j] : i < p < j}, 1 ≤ i < j-1 ≤ N-1 .

j j+1

Theorem : The length of the best 2-bridge move is min {g(i,j): 1 ≤ i < j-1 ≤ N-2 .

p+1 p

Question : How is this O(N 2 ) ?

There are and θ(N 2 ) values of f(p,j) g(i,j) time θ(N) to compute, and the average ones appears to take to compute.

i+1 i

2 1

N N-1 q+1 q

## The Functions to be Computed

• • f(p,j) = min {C[p,q] : j < q ≤ N}, 2 ≤ p < j ≤ N .

g(i,j) = C[i,j] + min {f[p,j] : i < p < j}, 1 ≤ i < j-1 ≤ N-1 .

j j+1

Observation : The only difference between f(p,j) and f(p,j-1) of C[p,j] is the inclusion in the minimization.

The only difference between g(i,j) and g(i-1,j) the inclusion of f(i,j) minimization.

is in the So the values of each can be computed in constant time per value.

p+1 p i+1 i

2 1

N N-1 q+1 q

## 4-Opt Conclusions

• • • • • Can be implemented to run in O(N 2 ) per iteration.

This is likely to be significantly slower than our 3-opt implementation.

k-opt for k > 4 will be even worse: – Intermediate solutions can involve a path plus multiple cycles.

– More possibilities for non-sequential moves that must be considered as special cases (triple-bridge moves, etc.), Improvement in tours over what 3-opt can produce may not be worth the effort (Shen Lin, private communication).

Alternative approach: a highly-pruned N-opt approach:

## The Lin-Kernighan Algorithm

### [Shen Lin & Brian Kernighan, Bell Labs, 1971]

• • Built on top of 3-opt, augmented to allow choices of t 5 that create a short cycle, so long as there are legal choices of t 6 , t 7 , and t 8 that pull things together again.

After each choice of t 6 (or t 8 in in the above special case), we invoke a deep-dive “LK-search,” assuming the gain criterion from the Partial Sum Theorem is met.

The source of what follows is [Lin & Kernighan, “An effective heuristic algorithm for the traveling salesman problem,” Operations Research 21 (1973), 498-516] and the original Lin-Kernighan FORTRAN code.

### LK Search

t 1 t 2i+1 t 2i+2 t 2i

Candidates for t 2i t 2i+1 are the cities on the neighbor list for such that the length of the resulting “one-tree” is less than the length of the shortest tour seen so far.

[Note that this is the same criterion as provided by the Partial Sum Theorem.] We abandon the choice if the edge to be deleted, tour edge). This limits the depth of the search to {t N 2i+1 .

,t 2i+2 } , was added to the tour earlier in this search (is not an original

t 2i+2

### LK Search

t 1 t 2i+1 t 2i t 2i+2

Candidates for t 2i t 2i+1 {t 1 ,t 2i+2 } to this path yields a new champion tour, such that the length of the resulting “one-tree” is less than the length of the shortest tour seen so far.

[Note that this is the same criterion as provided by the Partial Sum Theorem.] t 2i+1 .

The flips illustrated here are actually performed, and then undone as needed, eventually returning us to our starting point (or the improved tour, if found).

### Differences between Algorithm Tested and Lin-Kernighan, as described

• • • • • • • We use pre-computed neighbor lists ( k = 20 ). LK did not similarly restrict t 3 , t 5 , t 7 , although their code restricted LK search to the 5 nearest neighbors. We use don’t-look-bits and queue order. LK cycled through t 1 order.

candidates in input We use randomized Greedy starting tours, whereas LK used random starting tours.

Lin-Kernighan not only forbids the deletion of a previously-added tour edge. It also forbids addition of a previously deleted edge. We allow this latter possibility.

Lin-Kernighan also added a search for an improving double-bridge move (one that does not deleted added edges), used only whenever no further standard improving moves could be found. They used an exhaustive search, but even our O(N 2 ) algorithm is unlikely to scale well, so we omit this step.

Lin-Kernighan used Array tour representation, while we switch to 2-level trees for N > 1,000 .

Lin-Kernighan coded in FORTRAN and ran on memory-constrained machines. The largest instance they could test had 318 cities and was considered “big” in 1971. We used C and modern machines, and go considerably bigger.

## Show Movie

### Worst-Case Running Time per Iteration

(Time to find the next improving move or determine there is none.) • • • • • • • • • • O(N) choices for t 1 .

Up to 2 choices for t 2 .

In practice much smaller because of don’t-look-bits.

O(k) choices for t 3 , each potentially involving a between + a flip.

Up to 2 choices for t 4 .

O(k) choices for t 5 , each potentially involving a between + a flip.

Up to 2 choices for t 6 .

O(k) choices for t 7 , each potentially involving a between + a flip.

Up to 2 choices for t 8 .

O(N) levels of LK-search.

Typically much smaller.

O(k) choices for t 2i+1 , each potentially involving a flip.

2

4

### logN) using splay trees.

Further restricted by the locking of edges

## How Many Iterations?

• • • • • • In practice, typically θ(N) opt, and Lin-Kernighan.

on random Euclidean instances for all of 2-opt, 3 In theory, what?

Theorem [Englert, Röglin, & Vöcking, “Worst case and probabilistic analysis of the 2-opt algorithm for the TSP,” Algorithmica 68 (2014), 190-264] 1 , there exists a set of 16N points in the unit square for which 2-opt can take as many as 2 N+4 -22 steps.

: For any L p metric, 1 ≤ p ≤ ∞ , and any N ≥ One can get similar exponential lower bounds for k-opt, k > 2, if one considers instances not obeying the triangle inequality. [Chandra, Karloff, & Tovey, “New results on the old k-opt algorithm for the TSP,” SIAM J. Comput., 28 (1999), 1998–2029] .

Note: This is doubly worst-case – not only must one be unfortunate enough to get the bad instance, one must also be unfortunate enough to pick the bad sequence of moves. What about average case, at least on the instances?

For random Euclidean instances, one of these levels of worst-case can be removed.

### Iterations for Random Euclidean Instances

Theorem [Kern, 1989]: With high probability, maximum length of a sequence of improving 2-opt moves is O(N 16 ) .

Theorem [Chandra, Karloff, & Tovey, 1999]: The expected length of a maximum sequence of 2-opt moves is O(N 10 logN) and O(N 6 logN) for the Manhattan ( L 1 ) metric.

Theorem [Englert, Röglin, & Vöcking, 2014]: The expected length of a maximum sequence of 2-opt moves is O(N 4+1/3 logN) and O(N 3.5

logN) for the Manhattan metric.

(This latter result extends parametrically to higher dimensions and a variety of other point distributions.)

### Removing the Move Sequence from the Picture: PLS-completeness

[Johnson, Papadimitriou, & Yannakakis, “How easy is local search?” J. Comp. Syst. Sci. 37 (1988), 79-100] Definition : A local search problem L in PLS (polynomial-time local search) consists of a) A type T L ε { min , max }.

b) A polynomial-time recognizable set D L of instances.

c) For each instance x 𝜖 polynomial in |x| .

D L , a set F L (x) of solutions that is recognizable in time d) For each solution s 𝜖 F L (x) , 1) a non-negative integer cost c L (s,x) , and a 2) a subset N(s,x) ⊆ F L (x) , called the neighborhood of x.

e) Three polynomial-time algorithms: 1) 2) 3) A L , which, given x 𝜖 D L , produces a standard (starting) solution A L (x) 𝜖 F L (x) .

B L , which, given x 𝜖 D L and s 𝜖 F L (x) , computes c L (s,x) .

C L , which, given x 𝜖 D L and s 𝜖 F L better cost (e,g., c L (s’,x) < c L (s,x) solution exists, and hence s (x) , either returns a solution is locally optimal.

s’ 𝜖 N(s,x) with a if T L = min ), or reports truthfully that no such

in

### PLS:

• • • Is there an algorithm that, given an instance x 𝜖 solution s 𝜖 F L (x) in time polynomial in |x| ?

D L , finds a locally optimal True if all sequences of improving moves are polynomially bounded, as they would be, for instance, if all costs are polynomially bounded in |x| . Examples include – – – Vertex Cover Max Clique TSP when all edge lengths are less than some constant (as in the case of our random Euclidean instance generator, where all coordinates are in [0, 10 7 ] and distances are rounded to the nearest integer) Also may be true if there is a better heuristic for choosing the next move than the one given by C L (s,x) .

Or if there is some way of finding a locally optimal solution without using a local search algorithm at all.

– As when L is corresponds to the Simplex neighborhood for Linear Programming.

## PLS-Reductions

A reduction from PLS problem L time computable functions: to PLS problem K consists of two polynomial 1) 2) f g , which maps instances of , which maps pairs (x 𝜖 x 𝜖 D L , s 𝜖 D L to instances F K (f(x))) f(x) 𝜖 to solutions D K and g(x,s) 𝜖 F L (x) , such that for all then g(x,s) x 𝜖 D L , if s is a locally optimal solution for instance is locally optimal for L .

f(x) of K , This notion is transitive, and has all the properties needed to define the concept of PLS-completeness, and yield the following: Theorem : If we can find locally optimal solutions in polynomial time for any PLS-complete problem, then we can do so for all problems in PLS.

In addition, the running time of the “standard algorithm” for a PLS-complete problem is “typically” exponential, and the problem of determining the output of that algorithm is “typically” NP-hard.

(This is a property of many particular proofs of PLS-completeness, rather than a known consequence of the definitions).

### PLS-Completeness and the TSP

• • [Krentel, “Structure in locally optimal solutions,” FOCS 1989, IEEE Computer Society, pp. 216-221]: k-opt is PLS-complete for some k between 1,000 and 10,000.

[Papadimtriou, “The complexity of the Lin-Kernighan heuristic for the traveling salesman problem,” SIAM J. Comput. 7 (1992), 450-465]. Lin Kernighan is PLS-complete.

This is for the variant in which, for LK-search, instead of not allowing an added edge to be subsequently deleted, we instead do not allow an edge that has been deleted to be subsequently added back.

Note that under this variant, the LK-search can take as many as θ(N 2 ) steps, as compared to the O(N) bound for the other method.

And recall that

both

criteria are used in the original Lin-Kernighan paper.

### Computational Results, Random Euclidean Instances

2-Opt  3-Opt  LK 

N =

% Excess over HK 150 Mhz Secs % Excess 150 Mhz Secs % Excess 150 Mhz Secs

10 3

4.9

0.32

3.1

0.38

2.0

0.77

10 4

5.0

3.8

3.0

4.6

2.0

9.8

10 5

4.9

56.7

3.0

66.1

2.0

151

10 6

4.9

928 3.0

1054 2.0

2650 Time on 3.06 Ghz Intel Core i3 processor at N = 10 6 : 25.4 sec (2-opt), 29.5 sec (3-opt), 61.5 (LK) Of this, 23.5 sec was for Preprocessing (input reading + neighbor list construction + initial tour generation)

### Beating Lin-Kernighan, Part 1: Simulated Annealing

[Kirkpatrick, Gelatt, & Vecchi, “Optimization by simulated annealing,” Science 220 (13 May 1983), 671-680] [Černy, “A thermodynamical approach to the travelling salesman problem: An efficient simulation algorithm,” J. Optimization Theory and Appl. 45 (1985), 41-51] [Kirkpatrick, “Optimization by simulated annealing: Quantitative studies,” J. Stat. Physics 34 (1984), 976-986] Based on an analogy:

Physical System

State Energy Ground State Rapid Quenching Careful Annealing

Optimization Problem

Feasible Solution Cost Optimal Solution Local Optimization Simulated Annealing

## Theoretical Results

• • •

### Implementations Details

• • • • • Initial temperature so high that most moves are accepted.

Exponentials are evaluated by lookup from a pre-computed table. “Frozen” = No more moves being accepted.

“At Equilibrium” = Having tried a given fixed number of moves at the current temperature.

“Lower the temperature” = Multiply it by a fixed constant, say 0.95.

A generic simulated annealing implementation reflecting these principles was developed by DSJ, together with interns Cecilia Aragon, Lyle McGeoch, and Cathy Schevon. We consulted with Scott Kirkpatrick and so that our implementation would reflect his own. It is described in [Johnson, Aragon, McGeoch, & Schevon, “Optimization by simulated annealing: an experimental evaluation; Part I, Graph Partitioning,” Oper. Res 37 (1989), 865-892] and [Johnson, Aragon, McGeoch, & Schevon, “Optimization by simulated annealing: an experimental evaluation; Part II, Graph coloring and number partitioning,” Oper. Res 39 (1991), 378-406].

The implementation was adapted to the TSP with Lyle McGeoch. Results are described in [Johnson & McGeoch, “The traveling salesman problem: A case study in local optimization,” in Local Search in Combinatorial Optimization, Aarts & Lenstra (editors), John-Wiley and Sons, Ltd., 1997, pp. 215-310].

## TSP-Specific Details

Basic move: 2-opt (as done by Kirkpatrick).

Initial tour: Random tour.

Starting temperature: set so 97% of moves are accepted.

Temperature length = N(N-1), approximately twice the neighborhood size.

Temperature reduction factor = 0.95.

Averages over 10 runs.

Times are on a 150 Mhz processor “Excess” is percent over HK bound

### Algorithm Engineering for Simulated Annealing

• • • Neighborhood Pruning (first suggested by [Bonomi & Lutton, 1984]).

In our version, we only consider moves corresponding to an arbitrary city as t cities on 1 , one of its tour neighbors as t 2 ’s neighbor list, a total of t 2 , and t chosen from the 20 possibilities..

This is further restricted as follows. Let neighbor of c c = t 1 t 3 such that either d(c,c”) ≤ d(c,c’) d(c,c”) – d(c,c’) and let c’ be the tour that is farther away. We only consider candidates c’’ for or the probability of accepting an is greater than 0.01.

At each temperature, we consider being chosen as a candidate.

20 α N candidate moves ( α ≥ 1 a parameter), so that each potential move has a reasonably probability of • Low-temperature starts (suggested by Kirkpatrick).

Start with a Nearest Neighbor tour and a temperature yielding an initial acceptance rate of about 50%. For random Euclidean instances we follow the suggestion of Bonomi and Lutton to use L/√N , where L is the length of our “unit” square.

For more structured, instances, however, SA LK by a factor of 5 or more), SA 2 2 can beat LK on a time equalized basis.

In particular, for dsj1000 from TSPLIB (a clustered Euclidean instance which causes has 1.27% excess and time equalized LK gets 1.35%.

dsj1000

### More Algorithm Engineering: Using 3-opt Moves

• • Proposed by Kirkpatrick in 1984, but his implementation did not exploit what we now know about fast 3-opt, and so his move set only contained moves where the segment moved was of length 10 or less (a generalization of Or-opt).

In our version, we consider both 2- and 3-opt moves, choosing 2-opts for the first 10 α N candidates at each temperature, and 3-opts for the last 10 α N .

A random 3-opt move is chosen as follows. Choose t 1 N possibilities, t 2 the (pruned) list of up to 20 possibilities, t 4 randomly from the randomly from the two possibilities, t 3 randomly from randomly from the two possibilities, t 5 randomly from the 20 possibilities (with up to four retries if the initial choices are topologically infeasible), and then, if t 6 is not topologically determined, choose it randomly from the two possibilities. This yields a total of at most 3200N possibilities.

## Other Metaheuristics

• • • •

### Neural Nets [Hopfield & Tank, 1985]

– Hopeless

– Parallel exhaustive search in a bathtub – does not scale

### Tabu Search [Glover, 1986]

– Based on an idea implicit in LK-search - not competitive for TSP

### Genetic Algorithms [Holland, 1975]

– Initial ideas were not competitive, but…

## Genetic Algorithm Schema

1. Generate an initial population solutions.

P consisting of a random set of k 2. While not yet converged, create a new generation of P as follows.

1) Select k’ distinct one- or two-item subsets of P (the mating strategy) 2) For each one-element subset, perform a random mutation to obtain a new solution.

3) For each two-element subset, perform a randomized crossover operation to obtain a new solution that reflects both parents.

4) Let P’ be the set of solutions generated by (2) and (3).

5) Using a selection strategy, choose k these survivors.

survivors from P ∪ P’ and replace P by 3. Return the best solution in P .

Standard “convergence” strategy: Stop if you have run for j improving the best solution.

generations without Standard “selection” strategy: Choose the k best solutions.

## Sample Crossover Operation

• • • Pick a segment S from tour A.

Delete all the cities in S from tour B (keeping the remaining cities in the same order).

Insert segment S into what is left of tour B.

If you are skeptical that operations like this can lead to a competitive algorithm, you are right.

A new idea is needed: The “Hybrid Genetic Algorithm” [Brady, “Optimization strategies gleaned biological evolution,” Nature 317 (1985), 804-806].

### Hybrid Genetic Algorithm Schema

1.

Generate an initial population solutions.

P consisting of a random set of k 2. Apply a given local optimization algorithm A letting the resulting local optimal solution s’ to each solution s replace s in P .

in P , 3. While not yet converged, create a new generation of P as follows.

1) Select k’ distinct one- or two-item subsets of P (the mating strategy) 2) For each one-element subset, perform a random mutation to obtain a new solution.

3) For each two-element subset, perform a randomized crossover operation to obtain a new solution that reflects both parents.

4) Let P’ be the set of solutions generated by (2) and (3), after first applying algorithm A to each.

5) Using a selection strategy, choose k survivors.

survivors from P ∪ P’ and replace P by these 4. Return the best solution in P .

New

+ We can use Lin-Kernighan as the local optimization algorithm. Presumably we will get better tours than for basic Lin-Kernighan.

For a reasonably population size, we will need to perform lots of LK’s. Is this the most cost-effective way of investing all that extra time?

We still have the problem of the ad-hoc crossover to contend with.

In 1989, Martin, Otto, & Felten suggested a way to remove both drawbacks. (Their paper was published as [“Large-step Markov chains for the TSP incorporating local search heuristics,” Complex Systems 5 (1991), 299-326] .)

• • •

### For the mutation, perform a random double bridge 4-opt move.

Note that, although this was originally viewed in the context of genetic algorithms, there is little left that is genetic about it. The common name for this approach (using Lin-Kernighan) is now – – “Iterated Lin-Kernighan” or “Chained Lin-Kernighan.”

• • • • • Minimal extra programming needed to turn a local search heuristic into an iterated local search heuristic.

Preprocessing is amortized (although this is true of multiple independent-run LK as well).

After the first run of LK, subsequent ones are all much faster, since they start with all but 8 don’t-look-bits still on.

The damage done to the tour by just changing 4 edges still turns out to open up significant room for LK to find new improvements.

Surprisingly good performance.

### Random Euclidean Instances

.96

1704 3.06 Ghz Intel Core i3 processor

### Selected TSPLIB Instances

Except for instance fl3795, ILK(N) is always within 0.05% – 0.25% of optimal.

### TSPLIB Instance fl3795

ILK(N), 20 nearest neighbors: ILK(N), 20 quad neighbors: 4.33% excess (1080 seconds, 3.06 Ghz processor) 1.09% excess ( 205 seconds, 3.06 Ghz processor)* *Average of three runs, one of which found the optimal and a second found optimal + 1.

### Beating Iterated Lin-Kernighan

• Chained Lin-Kernighan [Applegate, Cook, & Rohe, “Chained Lin-Kernighan for large traveling salesman problems,” INFORMS J. Comput. 15 (2003), 82-92 ] Broader and deeper search before LK-searches, More selective in choices of double-bridge moves (mixture of greedier choices, random walks, etc.) • Hybrid Genetic Algorithm [Nguyen, Yoshihara, Yamamori, and Yasunaga, “Implementation of an effective hybrid GA for large-scale traveling salesman problems,” IEEE Trans. Syst., Man, Cybernetics B37 (2007), 92 99] Genetic algorithm with crossovers + mutations + Lin-Kernighan, and large running times for bigger instances.

Helsgaun’s algorithm [Helsgaun, “An effective implementation of the Lin Kernighan traveling salesman heuristic,” European J. Oper. Res. 126 (2000), 106-130. Source code available at http://www.dat.ruk.dk/~keld/ ] • Different restart mechanism, broader initial search, … Tour merging [Cook & Seymour, “Tour merging via branch-decomposition,” INFORMS J. Comput. 15 (2003), 233-248]

### Random Distance Matrices

N = 1000 3162 10,000 N = 1000 3162 10,000 % Excess of HK bound Seconds on a 500 Mhz DEC Alpha Note: Helsgaun’s algorithm is the only heuristic we’ve seen that is not fazed by these instances.

### Heuristics for the Asymmetric TSP (a quick tour)

• • • • Tour Construction: Asymmetric Nearest Neighbor Asymmetric Greedy Match & Patch (Compute minimum-cost cycle cover, repeatedly patch together the two largest cyles) Contract or Patch • • • • • • Repeated Assignment Zhang’s Algorithm (Truncated Branch & Bound) Local Optimization: Asymmetric 3-Opt Iterated 3-Opt Kanellakis-Papadimitriou (Mimics Lin-Kernighan, and also uses “best double-bridge move” as a component) Helsgaun (Convert to symmetric instance and apply Helsgaun’s algorithm.

## Contract-or-Patch

Once all 2-cycles are removed, proceed as in Match & Patch.

### Repeated Assignment

Delete all but one city from each cycle, find new matching, and repeat.

The union of all the matching edges is a connected graph in which every vertex has equal in- and out-degrees. So there is an Euler tour. If we assume the involves no more than half the previous number of vertices, we have: Δ Inequality, the Euler tour can be traversed using shortcuts at no extra cost. Each matching is no longer than an optimal tour, again by the triangle inequality. Since each matching Theorem [Frieze, Galbiati, & Maffioli (1982]: Assuming Δ Inequality, RA(I) ≤ log(N)OPT(I) .

### Heuristics for the Asymmetric TSP (a quick tour)

• • • • Tour Construction: Asymmetric Nearest Neighbor Asymmetric Greedy Match & Patch (Compute minimum-cost cycle cover, repeatedly patch together the two largest cyles) Contract or Patch • • • • • • Repeated Assignment Zhang’s Algorithm (Truncated Branch & Bound) Local Optimization: Asymmetric 3-Opt (Only 3-opt moves that cause no reversals) Iterated Asymmetric 3-Opt Kanellakis-Papadimitriou (Mimics Lin-Kernighan, and also uses “best double-bridge move” as a component) Helsgaun (Convert to symmetric instance and apply Helsgaun’s algorithm.

### Results for Instance Generators of

[Cirasella, Johnson, McGeoch, & Zhang, “The asymmetric traveling salesman problem: Algorithms, instance generators, and tests,” ALENEX 2001, Lecture Notes in Computer Science 2153, Springer, pp. 32-59]