Optimized TSP Algorithms

Question

I am interested in ways to improve or come up with algorithms that are able to solve the Travelling salesman problem for about n = 100 to 200 cities.

The wikipedia link I gave lists various optimizations, but it does so at a pretty high level, and I don't know how to go about actually implementing them in code.

There are industrial strength solvers out there, such as Concorde, but those are way too complex for what I want, and the classic solutions that flood the searches for TSP all present randomized algorithms or the classic backtracking or dynamic programming algorithms that only work for about 20 cities.

So, does anyone know how to implement a simple (by simple I mean that an implementation doesn't take more than 100-200 lines of code) TSP solver that works in reasonable time (a few seconds) for at least 100 cities? I am only interested in exact solutions.

You may assume that the input will be randomly generated, so I don't care for inputs that are aimed specifically at breaking a certain algorithm.

Solution 1

200 lines and no libraries is a tough constraint. The advanced solvers use branch and bound with the Held–Karp relaxation, and I'm not sure if even the most basic version of that would fit into 200 normal lines. Nevertheless, here's an outline.

Held Karp

One way to write TSP as an integer program is as follows (Dantzig, Fulkerson, Johnson). For all edges e, constant w_e denotes the length of edge e, and variable x_e is 1 if edge e is on the tour and 0 otherwise. For all subsets S of vertices, ∂(S) denotes the edges connecting a vertex in S with a vertex not in S.

minimize sum_{edges e} w_e x_e
subject to
1. for all vertices v, sum_{edges e in ∂({v})} x_e = 2
2. for all nonempty proper subsets S of vertices, sum_{edges e in ∂(S)} x_e ≥ 2
3. for all edges e in E, x_e in {0, 1}

Condition 1 ensures that the set of edges is a collection of tours. Condition 2 ensures that there's only one. (Otherwise, let S be the set of vertices visited by one of the tours.) The Held–Karp relaxation is obtained by making this change.

3. for all edges e in E, x_e in {0, 1}
3. for all edges e in E, 0 ≤ x_e ≤ 1

Held–Karp is a linear program but it has an exponential number of constraints. One way to solve it is to introduce Lagrange multipliers and then do subgradient optimization. That boils down to a loop that computes a minimum spanning tree and then updates some vectors, but the details are sort of involved. Besides "Held–Karp" and "subgradient (descent|optimization)", "1-tree" is another useful search term.

(A slower alternative is to write an LP solver and introduce subtour constraints as they are violated by previous optima. This means writing an LP solver and a min-cut procedure, which is also more code, but it might extend better to more exotic TSP constraints.)

Branch and bound

By "partial solution", I mean an partial assignment of variables to 0 or 1, where an edge assigned 1 is definitely in the tour, and an edge assigned 0 is definitely out. Evaluating Held–Karp with these side constraints gives a lower bound on the optimum tour that respects the decisions already made (an extension).

Branch and bound maintains a set of partial solutions, at least one of which extends to an optimal solution. The pseudocode for one variant, depth-first search with best-first backtracking is as follows.

let h be an empty minheap of partial solutions, ordered by Held–Karp value
let bestsolsofar = null
let cursol be the partial solution with no variables assigned
loop
    while cursol is not a complete solution and cursol's H–K value is at least as good as the value of bestsolsofar
        choose a branching variable v
        let sol0 be cursol union {v -> 0}
        let sol1 be cursol union {v -> 1}
        evaluate sol0 and sol1
        let cursol be the better of the two; put the other in h
    end while
    if cursol is better than bestsolsofar then
        let bestsolsofar = cursol
        delete all heap nodes worse than cursol
    end if
    if h is empty then stop; we've found the optimal solution
    pop the minimum element of h and store it in cursol
end loop

The idea of branch and bound is that there's a search tree of partial solutions. The point of solving Held–Karp is that the value of the LP is at most the length OPT of the optimal tour but also conjectured to be at least 3/4 OPT (in practice, usually closer to OPT).

The one detail in the pseudocode I've left out is how to choose the branching variable. The goal is usually to make the "hard" decisions first, so fixing a variable whose value is already near 0 or 1 is probably not wise. One option is to choose the closest to 0.5, but there are many, many others.

EDIT

Java implementation. 198 nonblank, noncomment lines. I forgot that 1-trees don't work with assigning variables to 1, so I branch by finding a vertex whose 1-tree has degree >2 and delete each edge in turn. This program accepts TSPLIB instances in EUC_2D format, e.g., eil51.tsp and eil76.tsp and eil101.tsp and lin105.tsp from http://www2.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp/.

// simple exact TSP solver based on branch-and-bound/Held--Karp
import java.io.*;
import java.util.*;
import java.util.regex.*;

public class TSP {
  // number of cities
  private int n;
  // city locations
  private double[] x;
  private double[] y;
  // cost matrix
  private double[][] cost;
  // matrix of adjusted costs
  private double[][] costWithPi;
  Node bestNode = new Node();

  public static void main(String[] args) throws IOException {
    // read the input in TSPLIB format
    // assume TYPE: TSP, EDGE_WEIGHT_TYPE: EUC_2D
    // no error checking
    TSP tsp = new TSP();
    tsp.readInput(new InputStreamReader(System.in));
    tsp.solve();
  }

  public void readInput(Reader r) throws IOException {
    BufferedReader in = new BufferedReader(r);
    Pattern specification = Pattern.compile("\\s*([A-Z_]+)\\s*(:\\s*([0-9]+))?\\s*");
    Pattern data = Pattern.compile("\\s*([0-9]+)\\s+([-+.0-9Ee]+)\\s+([-+.0-9Ee]+)\\s*");
    String line;
    while ((line = in.readLine()) != null) {
      Matcher m = specification.matcher(line);
      if (!m.matches()) continue;
      String keyword = m.group(1);
      if (keyword.equals("DIMENSION")) {
        n = Integer.parseInt(m.group(3));
        cost = new double[n][n];
      } else if (keyword.equals("NODE_COORD_SECTION")) {
        x = new double[n];
        y = new double[n];
        for (int k = 0; k < n; k++) {
          line = in.readLine();
          m = data.matcher(line);
          m.matches();
          int i = Integer.parseInt(m.group(1)) - 1;
          x[i] = Double.parseDouble(m.group(2));
          y[i] = Double.parseDouble(m.group(3));
        }
        // TSPLIB distances are rounded to the nearest integer to avoid the sum of square roots problem
        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
            double dx = x[i] - x[j];
            double dy = y[i] - y[j];
            cost[i][j] = Math.rint(Math.sqrt(dx * dx + dy * dy));
          }
        }
      }
    }
  }

  public void solve() {
    bestNode.lowerBound = Double.MAX_VALUE;
    Node currentNode = new Node();
    currentNode.excluded = new boolean[n][n];
    costWithPi = new double[n][n];
    computeHeldKarp(currentNode);
    PriorityQueue<Node> pq = new PriorityQueue<Node>(11, new NodeComparator());
    do {
      do {
        boolean isTour = true;
        int i = -1;
        for (int j = 0; j < n; j++) {
          if (currentNode.degree[j] > 2 && (i < 0 || currentNode.degree[j] < currentNode.degree[i])) i = j;
        }
        if (i < 0) {
          if (currentNode.lowerBound < bestNode.lowerBound) {
            bestNode = currentNode;
            System.err.printf("%.0f", bestNode.lowerBound);
          }
          break;
        }
        System.err.printf(".");
        PriorityQueue<Node> children = new PriorityQueue<Node>(11, new NodeComparator());
        children.add(exclude(currentNode, i, currentNode.parent[i]));
        for (int j = 0; j < n; j++) {
          if (currentNode.parent[j] == i) children.add(exclude(currentNode, i, j));
        }
        currentNode = children.poll();
        pq.addAll(children);
      } while (currentNode.lowerBound < bestNode.lowerBound);
      System.err.printf("%n");
      currentNode = pq.poll();
    } while (currentNode != null && currentNode.lowerBound < bestNode.lowerBound);
    // output suitable for gnuplot
    // set style data vector
    System.out.printf("# %.0f%n", bestNode.lowerBound);
    int j = 0;
    do {
      int i = bestNode.parent[j];
      System.out.printf("%f\t%f\t%f\t%f%n", x[j], y[j], x[i] - x[j], y[i] - y[j]);
      j = i;
    } while (j != 0);
  }

  private Node exclude(Node node, int i, int j) {
    Node child = new Node();
    child.excluded = node.excluded.clone();
    child.excluded[i] = node.excluded[i].clone();
    child.excluded[j] = node.excluded[j].clone();
    child.excluded[i][j] = true;
    child.excluded[j][i] = true;
    computeHeldKarp(child);
    return child;
  }

  private void computeHeldKarp(Node node) {
    node.pi = new double[n];
    node.lowerBound = Double.MIN_VALUE;
    node.degree = new int[n];
    node.parent = new int[n];
    double lambda = 0.1;
    while (lambda > 1e-06) {
      double previousLowerBound = node.lowerBound;
      computeOneTree(node);
      if (!(node.lowerBound < bestNode.lowerBound)) return;
      if (!(node.lowerBound < previousLowerBound)) lambda *= 0.9;
      int denom = 0;
      for (int i = 1; i < n; i++) {
        int d = node.degree[i] - 2;
        denom += d * d;
      }
      if (denom == 0) return;
      double t = lambda * node.lowerBound / denom;
      for (int i = 1; i < n; i++) node.pi[i] += t * (node.degree[i] - 2);
    }
  }

  private void computeOneTree(Node node) {
    // compute adjusted costs
    node.lowerBound = 0.0;
    Arrays.fill(node.degree, 0);
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) costWithPi[i][j] = node.excluded[i][j] ? Double.MAX_VALUE : cost[i][j] + node.pi[i] + node.pi[j];
    }
    int firstNeighbor;
    int secondNeighbor;
    // find the two cheapest edges from 0
    if (costWithPi[0][2] < costWithPi[0][1]) {
      firstNeighbor = 2;
      secondNeighbor = 1;
    } else {
      firstNeighbor = 1;
      secondNeighbor = 2;
    }
    for (int j = 3; j < n; j++) {
      if (costWithPi[0][j] < costWithPi[0][secondNeighbor]) {
        if (costWithPi[0][j] < costWithPi[0][firstNeighbor]) {
          secondNeighbor = firstNeighbor;
          firstNeighbor = j;
        } else {
          secondNeighbor = j;
        }
      }
    }
    addEdge(node, 0, firstNeighbor);
    Arrays.fill(node.parent, firstNeighbor);
    node.parent[firstNeighbor] = 0;
    // compute the minimum spanning tree on nodes 1..n-1
    double[] minCost = costWithPi[firstNeighbor].clone();
    for (int k = 2; k < n; k++) {
      int i;
      for (i = 1; i < n; i++) {
        if (node.degree[i] == 0) break;
      }
      for (int j = i + 1; j < n; j++) {
        if (node.degree[j] == 0 && minCost[j] < minCost[i]) i = j;
      }
      addEdge(node, node.parent[i], i);
      for (int j = 1; j < n; j++) {
        if (node.degree[j] == 0 && costWithPi[i][j] < minCost[j]) {
          minCost[j] = costWithPi[i][j];
          node.parent[j] = i;
        }
      }
    }
    addEdge(node, 0, secondNeighbor);
    node.parent[0] = secondNeighbor;
    node.lowerBound = Math.rint(node.lowerBound);
  }

  private void addEdge(Node node, int i, int j) {
    double q = node.lowerBound;
    node.lowerBound += costWithPi[i][j];
    node.degree[i]++;
    node.degree[j]++;
  }
}

class Node {
  public boolean[][] excluded;
  // Held--Karp solution
  public double[] pi;
  public double lowerBound;
  public int[] degree;
  public int[] parent;
}

class NodeComparator implements Comparator<Node> {
  public int compare(Node a, Node b) {
    return Double.compare(a.lowerBound, b.lowerBound);
  }
}

Solution 2

If your graph satisfy the triangle inequality and you want a guarantee of 3/2 within the optimum I suggest the christofides algorithm. I've wrote an implementation in php at phpclasses.org.

Solution 3

As of 2013, It is possible to solve for 100 cities using only the exact formulation in Cplex. Add degree equations for each vertex, but include subtour-avoiding constraints only as they appear. Most of them are not necessary. Cplex has an example on this.

You should be able to solve for 100 cities. You will have to iterate every time a new subtour is found. I ran an example here and in a couple of minutes and 100 iterations later I got my results.

Solution 4

I took Held-Karp algorithm from concorde library and 25 cities are solved in 0.15 seconds. This performance is perfectly good for me! You can extract the code (writen in ANSI C) of held-karp from concorde library: http://www.math.uwaterloo.ca/tsp/concorde/downloads/downloads.htm. If the download has the extension gz, it should be tgz. You might need to rename it. Then you should make little ajustments to port in in VC++. First take the file heldkarp h and c (rename it cpp) and other about 5 files, make adjustments and it should work calling CCheldkarp_small(...) with edgelen: euclid_ceiling_edgelen.

View more solutions

Comments

[Micromega]

Exact solutions and 100 cities? Are you mad?

[IVlad]

@Jitamaro - I would say that yes, I am mad, but I don't see how that's relevant here.

[Micromega]

I know some good heuristics but you wrote EXACT SOLUTIONS. You should correct your question.

[IVlad]

@Jitamaro - no, I mean Exact Solutions. Heuristics can also mean that you do not explore paths that will definitely not lead to a solution, not necessarily that you will settle for an approximation. See the "exact algorithms" section on wikipedia.

[IVlad]

I know it's NP-hard. That alone doesn't mean it won't finish for 100. There are various heuristics you can apply so that you eliminate a lot of possibilities that CANNOT generate a solution. Please read the wikipedia article for summary descriptions of those methods. If you don't trust wikipedia, there are also books that mention them, and existing programs that work fast for 100 cities. I am just asking how. NP-hard just means that no polynomial algorithm exists. It doesn't mean there's no exact algorithm that finishes in 1 second for n = 100.

[Karoly Horvath]

I hope you're talking about this node only: en.wikipedia.org/wiki/… . the rest is just heuristics which get "close" to an optimal solution.

[IVlad]

yes, only about that node, but I did read something about how the k-opt heuristics can yield exact solutions once. I don't have it anymore and I don't know if it was right or not, just saying.

[Karoly Horvath]

can yield != will yield. so it's useless if you are looking for an optimal solution.

[IVlad]

By can yield I meant that it might be possible to adapt it so it will always give an optimal solution.

[IVlad]

I would like a general solution for any type of graph, but I didn't know about that algorithm, so thanks. I think you should leave this up.

[IVlad]

I don't see why it would take that long for just 100 nodes. The Wikipedia article states solvers that have solved instances with tens of thousands of cities - on supercomputers and in 24 hours, sure, but reducing their input to 100 cities makes it seem like they would do it pretty fast on a single processor. What about the techniques that Wikipedia mentions for 100 - 200 instances? Do you have more details on those? I still see no reason why this shouldn't be possible. Everyone keeps saying it's because it's NP hard, but that's really just a theoretical result and doesn't prove anything.

[IVlad]

As for my exact figures - it's not like I expect someone to post an algorithm with exactly 200 lines of code that does it. They're just guidelines. I'm only interested in more details about the methods Wikipedia lists under exact algorithms, or other methods that achieve similar results.

[amit]

it's NP-Hard, not "as far as we know", it just is. if P=NP, it is still NP-Hard, because there is a reduction from every NP problem to TSP [definition of NP-Hard].

[Karoly Horvath]

@amit: I said it's NP-Hard. I think you missed the period between the sentences.

[Rob Neuhaus]

You probably won't be able to solve all n=100 problems exactly, but people have exactly solved a real instance of size n=33,810 (see wikipedia).

[flolo]

@rrenaud: yes, but who knows which problem specific facts they have used. I checked the references in the article and the Concorde numbers (btw. its from the same author whose book I recommend) and he can for 200 city examples indeed make it in 10s of seconds solve (I dont checked the examples, maybe there were also some problem specific knowledge used for optimization), but Concorde is not 200 lines of code, but 200 thousand (author states 150k). So I still say: 200 lines, exact, 200 nodes, few seconds - you cant get all 4 at the same time.

[IVlad]

+1, this is the kind of answer I was expecting. I will try to digest it for a few days and give more people the change to contribute, then I will most likely approve it.

[thomasfedb]

Nice answer, and your first too!

[Alexander]

Can you please explain what the input and output files format means. The example of input file I took from here github.com/pdrozdowski/TSPLib.Net/blob/master/TSPLIB95/tsp/…