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.

Answer 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);
  }
}