Report

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 Outline 1. k-Opt, k > 3 2. Lin-Kernighan 3. and Beyond…. 4. A quick tour of heuristics for the Asymmetric TSP 4-Opt • Seems to be relatively straightforward to generalize from 3-opt. – Just need to now consider possibilities for t7 and t8. – 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 t1 t3 t4 t2 t4 An alternative View: Maintain the path t1 t3 t2 t4 An alternative View: Maintain the path, at least eventually… t1 t4 t3 t5 t6 t2 An alternative View: Maintain the path, at least eventually… t1 t4 t35 t53 t26 t62 Double Bridge Move Cannot be produced by any sequential move. Best one can be found in O(N2) time. Best Double Bridge Move in Time O(N2) [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+1 p Smallest index of the first endpoint of a deleted edge q qp+1 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 Theorem: The length of the best 2-bridge move is min {g(i,j): 1 ≤ i < j-1 ≤ N-2. j+1 p+1 p Question: How is this O(N2)? There are values of f(p,j) and g(i,j) to compute, and the average ones appears to take time θ(N) to compute. q θ(N2) q+1 i+1 i 2 1 N N-1 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 Observation: The only difference between f(p,j) and f(p,j-1) is the inclusion of C[p,j] in the minimization. p+1 p The only difference between g(i,j) and g(i-1,j) is the inclusion of f(i,j) in the minimization. So the values of each can be computed in constant time per value. j+1 q q+1 i+1 i 2 1 N N-1 4-Opt Conclusions • Can be implemented to run in O(N2) 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 The Lin-Kernighan Algorithm [Shen Lin & Brian Kernighan, Bell Labs, 1971] • Built on top of 3-opt, augmented to allow choices of t5 that create a short cycle, so long as there are legal choices of t6, t7, and t8 that pull things together again. • After each choice of t6 (or t8 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 t1 t2i+1 t2i+2 t2i Candidates for t2i+1 are the cities on the neighbor list for t2i 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, {t2i+1,t2i+2}, was added to the tour earlier in this search (is not an original tour edge). This limits the depth of the search to N. t2i+2 LK Search t1 t2i+1 t2i t2i+2 If adding the Candidates foredge t2i+1 {t are the} cities to thison path theyields neighbor a new listchampion for tour, 1,t2i+2 t2i such save that that tour. the length of the resulting “one-tree” is less than the length of the shortest tour seen so far. As thethat starting thecriterion next level the search, take the [Note this path is thefor same asof provided by the shortest path generated by any of the the valid choices for t2i+1. Partial Sum Theorem.] We the choice if are edgeactually to be deleted, {t2i+1and ,t2i+2then }, undone The abandon flips illustrated here performed, was added to the tourreturning earlier (isusnot tour as needed, eventually to an ouroriginal starting point (or the edge). This limits the depth of the search to N. 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 t3, t5, t7, although their code restricted LK search to the 5 nearest neighbors. • We use don’t-look-bits and queue order. LK cycled through t1 candidates in input order. • 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(N2) 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 (Time to find the next improving move or determine there is none.) • O(N) choices for t1. • Up to 2 choices for t2. • O(k) choices for t3, each potentially involving a between + a flip. • Up to 2 choices for t4. • O(k) choices for t5, each potentially involving a between + a flip. • Up to 2 choices for t6. • O(k) choices for t7, each potentially involving a between + a flip. • Up to 2 choices for t8. • O(N) levels of LK-search. • O(k) choices for t2i+1, each potentially involving a flip. In practice much smaller because of don’t-look-bits. Typically much smaller. Total = O(N2k4logN) using splay trees. Further restricted by the locking of edges In practice much smaller, since tour edges usually go to near neighbors Worst-Case Running Time per Iteration How Many Iterations? • In practice, typically θ(N) on random Euclidean instances for all of 2-opt, 3opt, and Lin-Kernighan. • 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]: For any Lp metric, 1 ≤ p ≤ ∞, and any N ≥ 1, there exists a set of 16N points in the unit square for which 2-opt can take as many as 2N+4 -22 steps. • 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(N16). Theorem [Chandra, Karloff, & Tovey, 1999]: The expected length of a maximum sequence of 2-opt moves is O(N10logN) and O(N6logN) for the Manhattan (L1) metric. Theorem [Englert, Röglin, & Vöcking, 2014]: The expected length of a maximum sequence of 2-opt moves is O(N4+1/3logN) and O(N3.5logN) 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 TL ε {min, max}. b) A polynomial-time recognizable set DL of instances. c) For each instance x DL, a set FL(x) of solutions that is recognizable in time polynomial in |x|. d) For each solution s FL(x), 1) a non-negative integer cost cL(s,x), and a 2) a subset N(s,x) ⊆ FL(x), called the neighborhood of x. e) Three polynomial-time algorithms: 1) AL, which, given x DL, produces a standard (starting) solution AL(x) FL(x). 2) BL, which, given x DL and s FL(x), computes cL(s,x). 3) CL, which, given x DL and s FL(x), either returns a solution s’ N(s,x) with a better cost (e,g., cL(s’,x) < cL(s,x) if TL = min), or reports truthfully that no such solution exists, and hence s is locally optimal. Key Computational Question About a Problem L in PLS: Is there an algorithm that, given an instance x DL, finds a locally optimal solution s FL(x) in time polynomial in |x|? • 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, 107] 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 CL(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 to PLS problem K consists of two polynomialtime computable functions: 1) f, which maps instances of x DL to instances f(x) DK and 2) g, which maps pairs (x DL, s FK(f(x))) to solutions g(x,s) FL(x), such that for all x DL, if s is a locally optimal solution for instance f(x) of K, then g(x,s) is locally optimal for L. 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]. LinKernighan 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 θ(N2) 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 N= 103 104 105 106 2-Opt [20] % Excess over HK 4.9 5.0 4.9 4.9 0.32 3.8 56.7 928 3.1 3.0 3.0 3.0 0.38 4.6 66.1 1054 2.0 2.0 2.0 2.0 0.77 9.8 151 2650 150 Mhz Secs 3-Opt [20] % Excess 150 Mhz Secs LK [20] % Excess 150 Mhz Secs Time on 3.06 Ghz Intel Core i3 processor at N = 106: 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 Optimization Problem State Feasible Solution Energy Cost Ground State Optimal Solution Rapid Quenching Local Optimization Careful Annealing Simulated Annealing Theoretical Results General Theorem [Proved by many]: • If you run long enough, and cool slowly enough • (say letting the temperature be C/log(n), where C is a constant and n is the number of moves tested so far), • then, with high probability, you will converge to an optimal solution. General Drawback: For this to work, “long enough” will exceed the time to perform exhaustive search… 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 t1, one of its tour neighbors as t2, and t3 chosen from the 20 cities on t2’s neighbor list, a total of 20N possibilities.. • This is further restricted as follows. Let c = t1 and let c’ be the tour neighbor of c that is farther away. We only consider candidates c’’ for t3 such that either d(c,c”) ≤ d(c,c’) or the probability of accepting an uphill move of size d(c,c”) – d(c,c’) is greater than 0.01. • At each temperature, we consider 20αN candidate moves (α ≥ 1 a parameter), so that each potential move has a reasonably probability of being chosen as a candidate. 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, SA2 can beat LK on a time equalized basis. In particular, for dsj1000 from TSPLIB (a clustered Euclidean instance which causes LK by a factor of 5 or more), SA2 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 t1 randomly from the N possibilities, t2 randomly from the two possibilities, t3 randomly from the (pruned) list of up to 20 possibilities, t4 randomly from the two possibilities, t5 randomly from the 20 possibilities (with up to four retries if the initial choices are topologically infeasible), and then, if t6 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 • DNA computing [Adleman, 1994] – 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 P consisting of a random set of k solutions. 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 survivors from P∪P’ and replace P by these survivors. 3. Return the best solution in P. Standard “convergence” strategy: Stop if you have run for j generations without improving the best solution. Standard “selection” strategy: Choose the k best solutions. Sample Crossover Operation a b u h i q w o p v n m s z x l k e r y d f g j c q w e r t y u o p a d s f d g f h g j h c j v k b l n z m x s c z v x b l n k m • 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 P consisting of a random set of k solutions. 2. Apply a given local optimization algorithm A to each solution s in P, letting the resulting local optimal solution s’ replace s 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 from P∪P’ and replace P by these survivors. 4. Return the best solution in P. New Advantages and Disadvantages + 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].) Martin, Otto, & Felten’s Approach • Set population size to 1. • Only do mutations (no matings). • For the mutation, perform a random doublebridge 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.” Advantages of Iterated 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 multipleindependent-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 Random Euclidean Instances 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: 4.33% excess (1080 seconds, 3.06 Ghz processor) ILK(N), 20 quad neighbors: 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. Random Distance Matrices 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), 9299] • Genetic algorithm with crossovers + mutations + Lin-Kernighan, and large running times for bigger instances. Helsgaun’s algorithm [Helsgaun, “An effective implementation of the LinKernighan 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 % Excess of HK bound N = 1000 3162 10,000 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. Matching-and-Patch 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 Δ-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 involves no more than half the previous number of vertices, we have: 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] The Coin Collection Problem The Stacker-Crane Problem The No-Wait Shop Scheduling Problem The Disk Scheduling Problem The Shortest Common Superstring Problem The Drilling on a Tilted Table Problem Next Time Optimization Algorithms