TSP Approximation Algorithms, in depth: the double-tree bound, insertion heuristics, the Held-Karp 1-tree bound and measured gaps

Metric instances and why the triangle inequality matters

Every guarantee on this page assumes a metric instance: distances are symmetric, non-negative and obey the triangle inequality d(a, c) ≤ d(a, b) + d(b, c). Euclidean distances, shortest-path road distances and great-circle distances are metric. Raw travel-time tables often are not: one-way streets break symmetry, and time-dependent or toll-adjusted costs can break the triangle inequality.

The assumption is load-bearing. Without the triangle inequality there is no polynomial-time algorithm with any constant ratio unless P = NP (Sahni and Gonzalez, 1976). The reduction is short: give the edges of a graph cost 1 and the non-edges a huge cost. A tour of cost n exists exactly when the graph has a Hamiltonian cycle, so an approximation with any fixed ratio would decide an NP-complete problem. Even for metric instances, Karpinski, Lampis and Schmied showed in 2013 that approximating within a factor below 123/122 is NP-hard, so unless P = NP no polynomial-time approximation scheme exists for general metrics.

If your matrix is not metric, take its metric closure first: replace d(a, b) by the shortest-path distance from a to b (Floyd-Warshall, or Dijkstra from each city). A tour in the closure maps back to a walk in the original graph with the same cost, possibly passing through cities more than once, which is what a vehicle does anyway.

The tree lower bound and the double-tree 2-approximation

Delete any edge from an optimal tour and you get a Hamiltonian path, which is a spanning tree. So the minimum spanning tree weight satisfies MST ≤ OPT. That one inequality powers the simplest guarantee:

1. Build a minimum spanning tree T (Prim on a dense matrix in O(n²)).
2. Walk T depth first. The walk crosses every tree edge twice, so it costs 2·w(T).
3. Shortcut: list cities in the order they are first visited (the preorder) and jump straight from each to the next. By the triangle inequality every jump costs at most the stretch of walk it replaces.

The tour therefore costs at most 2·w(T) ≤ 2·OPT. The figure shows it on eight cities: the tree on the left, the preorder tour on the right, with dashed jumps where the walk would have doubled back. The ratio of 2 is tight for this method: there are metric instances where the shortcut tour approaches twice the optimum, and on random instances it is noticeably worse than other cheap methods (measured below).

Construction heuristics and what is proven about them

Rosenkrantz, Stearns and Lewis (1977) analysed the construction heuristics that everyone still uses. Their results, for metric instances, are the ones worth carrying around:

Nearest insertion earns its 2 by imitating Prim. The city it adds next is exactly the one Prim's algorithm would attach next, and inserting it between neighbours i and j on the tour costs at most twice the attaching edge by the triangle inequality. Summing over all insertions gives at most 2·MST. Nearest neighbour has no such link to the tree, and there are instances where its ratio grows logarithmically. Farthest insertion is the odd one out: no constant bound is known for it, yet on random instances it routinely builds the best initial tour of the cheap methods. Worst-case ratios and typical behaviour rank the methods differently, which is why you measure.

An implementation

The functions below work on a dense distance matrix D. They are deliberately plain Python so the logic is visible; a production version would vectorise the inner minimums or use neighbour lists.

import math

def tour_len(D, t):
    return sum(D[t[i]][t[(i + 1) % len(t)]] for i in range(len(t)))

def mst_prim(D, skip=None):
    verts = [v for v in range(len(D)) if v != skip]
    root = verts[0]
    best = {v: (D[root][v], root) for v in verts if v != root}
    parent = {root: None}
    while best:
        v = min(best, key=lambda u: best[u][0])
        parent[v] = best.pop(v)[1]
        for u in best:
            if D[v][u] < best[u][0]:
                best[u] = (D[v][u], v)
    return parent

def double_tree(D):
    parent = mst_prim(D)
    kids = {v: [] for v in parent}
    for v, par in parent.items():
        if par is not None:
            kids[par].append(v)
    tour, stack = [], [0]
    while stack:                          # preorder = Euler walk with shortcuts
        v = stack.pop()
        tour.append(v)
        stack.extend(reversed(kids[v]))
    return tour

def nearest_neighbor(D, start=0):
    tour, left = [start], set(range(len(D))) - {start}
    while left:
        v = min(left, key=lambda u: D[tour[-1]][u])
        tour.append(v)
        left.remove(v)
    return tour

def insertion(D, rule="nearest"):
    tour = [0]
    dmin = {v: D[0][v] for v in range(1, len(D))}   # distance to the tour
    while dmin:
        pick = min if rule == "nearest" else max
        k = pick(dmin, key=dmin.get)
        del dmin[k]
        if len(tour) == 1:
            tour.append(k)
        else:
            def added(i):
                a, b = tour[i], tour[(i + 1) % len(tour)]
                return D[a][k] + D[k][b] - D[a][b]
            tour.insert(min(range(len(tour)), key=added) + 1, k)
        for v in dmin:
            dmin[v] = min(dmin[v], D[k][v])
    return tour

Insertion keeps a running distance from each outside city to the tour, so choosing the next city is O(n) and the whole construction is O(n²). For local improvement, run 2-opt (reverse a segment whenever two crossing edges can be uncrossed) as shown in the Christofides article, and add Or-opt (move a run of one to three cities elsewhere) if you have time budget left.

Measuring the gap: the Held-Karp 1-tree bound

A 1-tree is a spanning tree on cities 1..n-1 plus the two cheapest edges from city 0. Every tour is a 1-tree in which every city has degree 2, so the minimum 1-tree is a lower bound. Held and Karp's trick is to add a penalty π(v) to every edge touching v. Every tour's cost rises by exactly 2·Σπ because each city has degree two, so the minimum 1-tree minus 2·Σπ is still a lower bound for every choice of π. Subgradient ascent then pushes π up on cities whose 1-tree degree exceeds 2 and down on leaves, so the tree is driven towards a tour shape:

def one_tree_bound(D, pi):
    n = len(D)
    W = [[D[i][j] + pi[i] + pi[j] for j in range(n)] for i in range(n)]
    parent = mst_prim(W, skip=0)
    deg, w = [0] * n, 0.0
    for v, par in parent.items():
        if par is not None:
            w += W[v][par]; deg[v] += 1; deg[par] += 1
    for v in sorted(range(1, n), key=lambda v: W[0][v])[:2]:
        w += W[0][v]; deg[0] += 1; deg[v] += 1
    return w - 2 * sum(pi), deg

def held_karp_bound(D, upper, iters=300):
    pi, best, lam = [0.0] * len(D), float("-inf"), 2.0
    for k in range(iters):
        lb, deg = one_tree_bound(D, pi)
        best = max(best, lb)
        g = [d - 2 for d in deg]
        norm = sum(x * x for x in g)
        if norm == 0:
            break                         # the 1-tree is a tour: optimal
        step = lam * (upper - lb) / norm  # Polyak-style step toward the best tour
        pi = [a + step * x for a, x in zip(pi, g)]
        if k % 20 == 19:
            lam /= 2
    return best

The maximum over π is the Held-Karp bound, equal to the value of the subtour elimination linear programme. It is usually close to the optimum, which is what makes it useful: the gap between your tour and this bound is an honest upper bound on how much you are leaving on the table. In our brute-force checks on 300 random instances of 4-8 cities the bound never exceeded the optimum, and was never more than 0.2% below it.

Worked run: 200 cities

The same code on 200 uniform random cities in a 1000 × 1000 square (seed 7). The Held-Karp bound came out at 10,819.6 and converged within 300 subgradient steps (1,000 steps moved it by 0.04). The MST weighed 9,554.8. Gaps are relative to the Held-Karp bound, so they slightly overstate the true gap to the optimum:

The double tree, with the best guarantee here, built the worst tour; farthest insertion, with no constant guarantee, built the best. Guarantees bound disasters; they do not rank typical behaviour. Local search closes most of the gap but not all, and without the bound you could not say any of this. Note also that the MST was about 12% below the 1-tree bound, so report the latter. On the small instances, worst ratios to the brute-force optimum were 1.40 for double tree, 1.36 for nearest neighbour, 1.22 for nearest insertion and 1.09 for farthest insertion, all inside their proven limits.

Asymmetric, Euclidean and non-metric variants

Asymmetric TSP is NP-hard to approximate below 75/74 (Karpinski, Lampis and Schmied), and the gap between that and the best algorithm is still wide. If your costs are one-way, check which guarantee you are actually relying on before quoting a ratio.

Operational guidance

• Validate the metric. Sample triples and test d(a, c) ≤ d(a, b) + d(b, c) plus symmetry. A single bad row from a geocoding error voids every guarantee on this page.
• Report the gap, not just the length. Log tour length, Held-Karp bound and their ratio per run. A sudden rise in the gap is a regression signal even when nobody knows the optimum.
• Use neighbour lists past a few thousand cities. The dense O(n²) matrix is 80 MB at n = 3,163 in float64. Restrict 2-opt and insertion to the k nearest neighbours (k around 10) from a k-d tree.
• Budget time explicitly. Construction is cheap; local search has a long tail. Stop on a wall-clock limit or when an improvement round gains less than a threshold, and keep the best tour so far.
• Use a real solver when the gap matters. LKH for near-optimal tours on large instances, Concorde for proven optima on instances it can handle.

Failure modes

• Non-metric data. Shortcuts can lengthen the tour; the 2 and 3/2 proofs silently stop applying.
• A bad step size. Too large a subgradient step oscillates; too small stalls. Track the best bound seen, not the last one, as the code does.
• Float equality in local search. Without a small tolerance, 2-opt can cycle on moves that change length by rounding noise.

Trade-offs

Choose by what you need to promise. If a contract or downstream system needs a hard worst-case bound, use Christofides or the double tree and report the bound. If you need good tours fast, farthest or nearest insertion plus 2-opt and Or-opt is a strong default. For near-optimal results on large instances, use a mature Lin-Kernighan implementation and spend your effort on data quality.

What to do next

1. Run the code on your own distance matrix after checking symmetry and the triangle inequality on random triples.
2. Compute the Held-Karp bound once and record the gap of your current production tours.
3. Swap your construction step for farthest insertion and measure the gap again, before and after local search.
4. Add neighbour lists and Or-opt, then re-measure time and gap at your real size.
5. Read the Christofides deep dive for the 3/2 proof, and the approximation algorithms overview for how the same lower-bound pattern works for other problems.
6. For instance sizes where you might afford exactness, compare against the Held-Karp dynamic programme up to about 20 cities, and read about approximation schemes to see why the Euclidean case admits one and the general metric case does not.