The Traveling Salesman Problem in Theory & Practice

1. The Traveling Salesman Problem in Theory & Practice Lecture 12: Optimal Tour Lengths for Random Euclidean Instances 15 April 2014 David S. Johnson dstiflerj@gmail.com http://davidsjohnson.net Seeley Mudd 523, Tuesdays and Fridays

2. Outline • The Beardwood-Halton-Hammersley Theorem • Empirical Approximations to Asymptopia • Properties of Optimal Tours • Student Presentation by Junde Huang on the missing details of the proof that the Gilmore-Gomory algorithm correctly runs in polynomial time.

3. Random Euclidean Instances(The Classic TSP Test Case) N = 10

4. Random Euclidean Instances(The Classic TSP Test Case) N = 10

5. Random Euclidean Instances(The Classic TSP Test Case) N = 100

6. Random Euclidean Instances(The Classic TSP Test Case) N = 1000

7. Random Euclidean Instances(The Classic TSP Test Case) N = 10000

8. N = 100,000

9. Random Euclidean Instances • Performance of heuristics and optimization algorithms on these instances are reasonably well-correlated with that for real-world geometric instances. • We can generate many samples, and very large instances. • Instances have convenient statistical properties.

10. Optimal Tour Lengths (Difference from Mean):One Million 100-City Instances -1e+07 -5e+06 0 +5e+06 Optimal Tour Lengths Appear to Be Normally Distributed

11. Optimal Tour Lengths (Difference from Mean):Ten Thousand 1000-City Instances -1e+07 -5e+06 0 +5e+06 With a standard deviation that appears to be independent of N

12. Well-Defined Tour Length Asymptotics • Expected optimal tour length for an N-city instance approaches CN for some constant C as N  . [Beardwood, Halton, and Hammersley, 1959] Key Open Question:What is the Value of C?

13. The Early History • 1959: BHH estimated C  .75, based on hand solutions for a 202-city and a 400-city instance. • 1977: Stein estimates C  .765, based on extensive simulations on 100-city instances. • Methodological Problems: • Not enough data • Probably not true optima for the data there is • Misjudges asymptopia

14. Stein: C = .765 BHH: C = .75 Figure from [Johnson, McGeoch, Rothberg, 1996]

15. What is the dependence on N ? • Expected distance to nearest neighbor proportional to 1/N, times n cities yields (N) • O(N) cities close to the boundary are missing some neighbors, for an added contribution proportional to (N)(1/N), or (1) • A constant number of cities are close to two boundaries (at the corners of the square), which may add an additional (1/N ) • This yields target function OPT/N = C + /N + /N

16. Asymptotic Upper Bound Estimates (Heuristic-Based Results Fitted to OPT/N = C + /N) • 1989:Ong & Huang estimate C ≤ .74, based on runs of 3-Opt. • 1994:Fiechter estimates C ≤ .73, based on runs of “parallel tabu search” • 1994: Lee & Choi estimate C ≤ .721, based on runs of “multicanonical annealing” • Still inaccurate, but converging? • Needed: A new idea.

17. New Idea (1995): Suppress the variance added by the “Boundary Effect”byusingToroidal Instances • Join left boundary of the unit square to the right boundary, top to the bottom. 

18. Toroidal Unit Ball

19. Toroidal Distance Computations

20. Toroidal Instance Advantages • No boundary effects. • Same asymptotic constant for E[OPT/N] as for planar instances [Jaillet, 1992] (although it is still only asymptotic). • Lower empirical variance for fixed N.

21. Toroidal Approaches 1996:Percus & Martin estimate C  .7120 ± .0002. 1996:Johnson, McGeoch, and Rothberg estimate C  .7124 ± .0002. 2004: Jacobsen, Read, and Saleur estimate C  .7119. Each coped with the difficulty of computing optima in a different way.

22. Percus-Martin(Go Small) • Toroidal Instances with N ≤ 100: • 250,000 samples, N = 12,13,14,15,16,17 (“Optimal” = best of 10 Lin-Kernighan runs) • 10,000 samples with N = 30 (“Optimal” = best of 5 runs of 10-step-Chained-LK) • 6,000 samples with N = 100 (“Optimal” = best of 20 runs of 10-step-Chained-LK) • Fit to OPT/N = (C + a/N + b/N2)/(1+1/(8N)) (Normalization by the expected distance to the kth nearest neighbor)

23. Jacobsen-Read-Saleur(Go Narrow) • Cities go uniformly on a 1 x 100,000 cylinder – that is, only join the top and bottom of the unit square and stretch the width by a factor of 100,000. • For W = 1,2,3,4,5,6, set N = 100,000W and generate 10 sample instances. • Optimize by using dynamic programming, where only those cities within distance k of the frontier (~kw cities) can have degree 0 or 1, k = 4,5,6,7,8. • Estimate true optimal for fixed W as k  . • Estimate unit square constant as W  . • With N ≥ 100,000, assume no need for asymptotics in N

24. Johnson-McGeoch-Rothberg(Go Held-Karp) • Observe that • the Held-Karp (subtour) bound also has an asymptotic constant, i.e., HK/n  CHK [Goemans, 1995] , and is easier to compute than the optimal. • (OPT-HK)/N has a substantially lower variance than either OPT or HK. • Estimate • CHK based on instances from N=100 to 316,228, using heuristics and Concorde-based error estimates • (C- CHK)based on instances with N = 100, 316, 1000, using Concorde for N ≤ 316 and Iterated Lin-Kernighan plus Concorde-based error estimates for N = 1000.

25. Concorde • Not only computes optimal tours – also can compute precise Held-Karp bounds. • Only a pre-release version was available in 1995 when Johnson-McGeoch-Rothberg was written. • Machines are much faster now, cycles are much cheaper, and Concorde is much better.

26. Running times (in seconds) for 10,000 Concorde runs on random 1000-city planar Euclidean instances (2.66 Ghz Intel Xeon processor in dual-processor PC, purchased late 2002. Range: 7.1 seconds to 38.3 hours

27. Running times (in seconds) for 1,000,000 Concorde runs on random 1000-city “Toroidal” Euclidean instances Range: 2.6 seconds to 6 hours

28. The New Data • Points chosen uniformly from a 10Mx10M grid • Solver: • Latest (2003) version of Concordewith a few bug fixes and adaptations for new metrics • Primary Random Number Generator: • RngStreampackage of Pierre L’Ecuyer. See • “AN OBJECT-ORIENTED RANDOM-NUMBER PACKAGE WITH MANY LONG STREAMS AND SUBSTREAMS,” Pierre L'ecuyer, Richard Simard, E. Jack Chen, W. David Kelton, Operations Research 50:6 (2002), 1073-1075

29. Toroidal Instances

30. Euclidean Instances

31. Standard Deviations N = 100 N = 1,000

32. 99% Confidence Intervals for OPT/Nfor Euclidean and Toroidal Instances

33. 99% Confidence Intervals for (OPT-HK)/Nfor Euclidean and Toroidal Instances

34. Gnuplot Least Squares fit for the Percus-Martin values of N --OPT/N = (C + a/N + b/N2)/(1+1/(8N)) C = .712234 ± .00017 versus originally claimed C = .7120 ± .0002

35. Least Squares fit for all data from [12,100] --OPT/N = (C + a/N + b/N2) C = .712333 ± .00006 versus C = .712234 ± .00017

36. Least Squares fit for toroidal data from [30,2000] -- OPT/N = (C + a/N + b/N2) C = .712401 ± .000005 versus C = .712333 ± .00006

37. What is the right function? Power Series in 1/N – (Suggested by Percus-Martin) Justification: Expected distance to the kth nearest neighbor is provably such a power series.

38. What is the right function? OPT/sqrt(N) = Power Series in 1/sqrt(N)) Justification: This is what we saw in the planar Euclidean case (although it was caused by boundaries).

39. What is the right function? OPT/sqrt(N) = = (1/sqrt(N) · (Power Series in 1/N)

40. What is the right function?

41. Effect of Data Range on Estimate[30,2000], [60,2000], [100,2000], [200,2000], [100,1000] 95% Confidence Intervals C + a/n.5 C + a/n.5 + b/n C + a/n.5 + b/n + c/n1.5 C + a/n + b/n2 + c/n3 C + a/n.5 + b/n.1.5 C + a/n.5 + b/n.1.5 + c/n2.5

42. The Winners? C + a/n + b/n2 + c/n3 C = .71240 ± .000005 C = .71240 ± .00002

43. Does the HK-based approach agree? Question

44. CHK = .707980 ± .000003 95% confidence interval derived using C + a/N + b/N2 functional form

45. C-CHK = .004419 ± .000002 95% confidence interval derived using C + a/N + b/N2 functional form

46. HK-Based Estimate C-CHK = .004419 ± .000002 + CHK = .707980 ± .000003 C = .712399 ± .000005 Versus (Conservative) Opt-Based Estimate C = .712400 ± .000020 Combined Estimate? C = .71240 ± .00001

47. “Explaining” The Expected Optimal Tour Length • The fraction of optimal tour edges that go to kth nearest neighbor seems to be going to a constant ak for each k.

48. Fraction of Optimal Tour Edges 1st Neighbor (44.6%) 2nd Neighbor (26.0%) 3rd Neighbor (13.6%) 4th Neighbor (7.1%) 5th Neighbor (3.9%) 6th Neighbor (2.1%) 7th Neighbor (1.2%) 8th Neighbor (0.66%) 9th Neighbor (0.37%) 10th Neighbor (0.21%) 11th Neighbor (0.12%) 19th Neighbor (.00014%) 20th Neighbor (.00008%)

49. “Explaining” The Expected Optimal Tour Length • The fraction of optimal tour edges that go to kth nearest neighbor seems to be going to a constant ak for each k. • If dk is the average distance from a city to its kth nearest neighbor, then dksqrt(N) also seems to be going to a constant for each k.

50. (√N)·(Average distance to kth Nearest Neighbor) 20th Neighbor 17th Neighbor 14th Neighbor 11th Neighbor 8th Neighbor 6th Neighbor 5th Neighbor 4th Neighbor 3rd Neighbor 2nd Neighbor 1st Neighbor