Travelling Salesman Problem: Exact and Approximate Solutions
The Travelling Salesman Problem asks a deceptively simple question: given a set of cities and distances between them, what's the shortest route that visits each city exactly once and returns to the...
Key Insights
- Exact TSP solutions using Held-Karp dynamic programming run in O(n²·2ⁿ) time, making them practical only for problems under 20-25 cities; beyond that, you need approximation algorithms.
- Combining greedy initialization with 2-opt local search typically achieves solutions within 5% of optimal in milliseconds, offering the best quality-to-speed ratio for most practical applications.
- Metaheuristics like simulated annealing shine when you need better-than-greedy solutions but can’t afford exact methods; proper temperature scheduling is more important than algorithm choice.
Introduction to TSP
The Travelling Salesman Problem asks a deceptively simple question: given a set of cities and distances between them, what’s the shortest route that visits each city exactly once and returns to the starting point? Despite its simple formulation, TSP is NP-hard—no known polynomial-time algorithm solves it exactly.
This matters beyond academic curiosity. Logistics companies route delivery trucks through TSP variants daily. Circuit designers minimize wire length on PCBs. DNA sequencing algorithms use TSP formulations to reconstruct genomes from fragments. When Amazon optimizes last-mile delivery or when your GPS calculates a multi-stop route, TSP algorithms are doing the heavy lifting.
The problem’s difficulty scales brutally. With 10 cities, there are 181,440 possible tours. With 20 cities, that explodes to over 60 quadrillion. This article walks through the spectrum of solutions—from exact algorithms that guarantee optimality to heuristics that trade precision for tractability.
Brute Force & Dynamic Programming (Exact Solutions)
The naive approach enumerates all (n-1)!/2 possible tours. For symmetric TSP (where distance A→B equals B→A), we fix the starting city and consider half the permutations. This O(n!) complexity becomes unusable beyond 12-13 cities.
The Held-Karp algorithm, developed independently by Bellman and by Held and Karp in 1962, uses dynamic programming to reduce complexity to O(n²·2ⁿ). The key insight: we don’t need to track the order of visited cities, only which cities we’ve visited and where we currently are.
from functools import lru_cache
from typing import List, Tuple
def held_karp(distances: List[List[float]]) -> Tuple[float, List[int]]:
"""
Solve TSP exactly using Held-Karp dynamic programming.
Returns (minimum_cost, optimal_path).
"""
n = len(distances)
if n <= 1:
return 0, list(range(n))
# dp[visited_mask][current_city] = (min_cost, previous_city)
@lru_cache(maxsize=None)
def dp(visited: int, current: int) -> Tuple[float, int]:
if visited == (1 << n) - 1: # All cities visited
return distances[current][0], 0
min_cost = float('inf')
best_next = -1
for next_city in range(n):
if visited & (1 << next_city): # Already visited
continue
new_visited = visited | (1 << next_city)
cost = distances[current][next_city] + dp(new_visited, next_city)[0]
if cost < min_cost:
min_cost = cost
best_next = next_city
return min_cost, best_next
# Reconstruct path
total_cost, _ = dp(1, 0) # Start at city 0
path = [0]
visited = 1
current = 0
while visited != (1 << n) - 1:
_, next_city = dp(visited, current)
path.append(next_city)
visited |= (1 << next_city)
current = next_city
return total_cost, path
Held-Karp works for problems up to about 20-25 cities on modern hardware. The space complexity of O(n·2ⁿ) often becomes the limiting factor before time does.
Branch and Bound
Branch and bound systematically explores the solution space while pruning branches that can’t improve on the best known solution. The algorithm maintains a lower bound on partial solutions; if this bound exceeds our current best complete solution, we abandon that branch.
For TSP, we can compute lower bounds using minimum spanning trees. Any tour is a spanning tree plus one edge, so MST cost provides a valid lower bound.
import heapq
from typing import List, Tuple, Set
def mst_lower_bound(distances: List[List[float]],
visited: Set[int],
current: int,
start: int = 0) -> float:
"""
Compute MST-based lower bound for remaining cities.
"""
remaining = [i for i in range(len(distances)) if i not in visited]
if not remaining:
return distances[current][start]
# Include edges from current to remaining and remaining to start
nodes = remaining + [current, start]
# Prim's algorithm for MST
if len(remaining) == 1:
return distances[current][remaining[0]] + distances[remaining[0]][start]
in_mst = {current}
edges = []
total = 0
for node in remaining:
heapq.heappush(edges, (distances[current][node], node))
while len(in_mst) < len(nodes) - 1 and edges:
cost, node = heapq.heappop(edges)
if node in in_mst:
continue
in_mst.add(node)
total += cost
for next_node in remaining:
if next_node not in in_mst:
heapq.heappush(edges, (distances[node][next_node], next_node))
# Add minimum edge back to start
min_return = min(distances[node][start] for node in remaining)
return total + min_return
def branch_and_bound_tsp(distances: List[List[float]]) -> Tuple[float, List[int]]:
"""
Solve TSP using branch and bound with MST lower bounds.
"""
n = len(distances)
best_cost = float('inf')
best_path = []
# Priority queue: (lower_bound, current_cost, current_city, path, visited)
initial_bound = mst_lower_bound(distances, {0}, 0)
pq = [(initial_bound, 0, 0, [0], {0})]
while pq:
bound, cost, current, path, visited = heapq.heappop(pq)
if bound >= best_cost:
continue
if len(path) == n:
total = cost + distances[current][0]
if total < best_cost:
best_cost = total
best_path = path[:]
continue
for next_city in range(n):
if next_city in visited:
continue
new_cost = cost + distances[current][next_city]
new_visited = visited | {next_city}
new_path = path + [next_city]
lower_bound = new_cost + mst_lower_bound(
distances, new_visited, next_city
)
if lower_bound < best_cost:
heapq.heappush(pq, (
lower_bound, new_cost, next_city, new_path, new_visited
))
return best_cost, best_path
Branch and bound’s performance varies wildly based on problem structure. Tight lower bounds and good initial solutions dramatically improve pruning. In practice, it can solve instances up to 40-60 cities when the geometry is favorable.
Greedy & Nearest Neighbor Heuristics
When exact solutions are infeasible, greedy heuristics provide fast approximations. The nearest neighbor algorithm starts at an arbitrary city and repeatedly visits the closest unvisited city.
import math
from typing import List, Tuple
def nearest_neighbor(distances: List[List[float]],
start: int = 0) -> Tuple[float, List[int]]:
"""
Construct tour using nearest neighbor heuristic.
Returns (tour_cost, path).
"""
n = len(distances)
visited = [False] * n
path = [start]
visited[start] = True
total_cost = 0
current = start
for _ in range(n - 1):
nearest = -1
nearest_dist = float('inf')
for city in range(n):
if not visited[city] and distances[current][city] < nearest_dist:
nearest = city
nearest_dist = distances[current][city]
path.append(nearest)
visited[nearest] = True
total_cost += nearest_dist
current = nearest
total_cost += distances[current][start]
return total_cost, path
def multi_start_nearest_neighbor(distances: List[List[float]]) -> Tuple[float, List[int]]:
"""
Run nearest neighbor from each starting city, return best result.
"""
best_cost = float('inf')
best_path = []
for start in range(len(distances)):
cost, path = nearest_neighbor(distances, start)
if cost < best_cost:
best_cost = cost
best_path = path
return best_cost, best_path
Nearest neighbor runs in O(n²) time and typically produces tours 20-25% longer than optimal. Starting from multiple cities and keeping the best result improves average performance with linear overhead.
2-Opt and Local Search Improvements
Greedy solutions often contain obvious inefficiencies—crossing edges that could be “uncrossed” to shorten the tour. The 2-opt algorithm systematically tests all edge swaps.
A 2-opt move removes two edges and reconnects the tour differently. If the tour visits cities A→B→…→C→D, we test replacing edges (A,B) and (C,D) with (A,C) and (B,D), reversing the segment between B and C.
from typing import List, Tuple
def calculate_tour_cost(distances: List[List[float]], path: List[int]) -> float:
"""Calculate total tour cost for a given path."""
cost = sum(distances[path[i]][path[i+1]] for i in range(len(path)-1))
return cost + distances[path[-1]][path[0]]
def two_opt(distances: List[List[float]],
path: List[int]) -> Tuple[float, List[int]]:
"""
Improve tour using 2-opt local search.
Continues until no improving swap exists.
"""
n = len(path)
improved = True
path = path[:]
while improved:
improved = False
for i in range(n - 1):
for j in range(i + 2, n):
# Skip adjacent edges (no valid swap)
if j == i + 1:
continue
# Current edges: (i, i+1) and (j, j+1 mod n)
# New edges: (i, j) and (i+1, j+1 mod n)
i_next = i + 1
j_next = (j + 1) % n
# Calculate change in tour length
current = (distances[path[i]][path[i_next]] +
distances[path[j]][path[j_next]])
new = (distances[path[i]][path[j]] +
distances[path[i_next]][path[j_next]])
if new < current - 1e-10: # Improvement found
# Reverse segment between i+1 and j
path[i_next:j+1] = reversed(path[i_next:j+1])
improved = True
return calculate_tour_cost(distances, path), path
def greedy_plus_two_opt(distances: List[List[float]]) -> Tuple[float, List[int]]:
"""
Combine nearest neighbor initialization with 2-opt refinement.
"""
_, initial_path = multi_start_nearest_neighbor(distances)
return two_opt(distances, initial_path)
2-opt typically improves greedy solutions by 5-15%. The combination of multi-start nearest neighbor plus 2-opt often achieves solutions within 2-5% of optimal in O(n³) worst-case time, though practical performance is usually much better.
Metaheuristics: Simulated Annealing & Genetic Algorithms
Local search gets trapped in local optima. Metaheuristics introduce controlled randomness to escape these traps. Simulated annealing, inspired by metallurgical annealing, occasionally accepts worse solutions with probability decreasing over time.
import random
import math
from typing import List, Tuple
def simulated_annealing_tsp(
distances: List[List[float]],
initial_temp: float = 10000,
cooling_rate: float = 0.9995,
min_temp: float = 1e-8,
iterations_per_temp: int = 100
) -> Tuple[float, List[int]]:
"""
Solve TSP using simulated annealing with 2-opt moves.
"""
n = len(distances)
# Start with nearest neighbor solution
_, current_path = nearest_neighbor(distances)
current_cost = calculate_tour_cost(distances, current_path)
best_path = current_path[:]
best_cost = current_cost
temp = initial_temp
while temp > min_temp:
for _ in range(iterations_per_temp):
# Generate neighbor using 2-opt move
i = random.randint(0, n - 2)
j = random.randint(i + 2, n - 1 if i > 0 else n)
# Calculate delta (change in cost)
i_next = i + 1
j_next = (j + 1) % n
delta = (
distances[current_path[i]][current_path[j]] +
distances[current_path[i_next]][current_path[j_next]] -
distances[current_path[i]][current_path[i_next]] -
distances[current_path[j]][current_path[j_next]]
)
# Accept or reject move
if delta < 0 or random.random() < math.exp(-delta / temp):
current_path[i_next:j+1] = reversed(current_path[i_next:j+1])
current_cost += delta
if current_cost < best_cost:
best_cost = current_cost
best_path = current_path[:]
temp *= cooling_rate
return best_cost, best_path
Temperature scheduling is critical. Start too hot and you waste iterations on random walks. Cool too fast and you get stuck in local optima. The cooling rate of 0.9995 works well for most instances; adjust based on problem size and available time.
Benchmarking & Choosing the Right Approach
Here’s a benchmark harness to compare algorithms on real data:
import time
from typing import Callable, Dict, List, Tuple
def generate_random_instance(n: int, seed: int = 42) -> List[List[float]]:
"""Generate random Euclidean TSP instance."""
random.seed(seed)
points = [(random.uniform(0, 1000), random.uniform(0, 1000)) for _ in range(n)]
distances = [[0.0] * n for _ in range(n)]
for i in range(n):
for j in range(n):
dx = points[i][0] - points[j][0]
dy = points[i][1] - points[j][1]
distances[i][j] = math.sqrt(dx*dx + dy*dy)
return distances
def benchmark_algorithms(sizes: List[int] = [10, 15, 20, 50, 100]) -> None:
"""
Compare all TSP algorithms across different problem sizes.
"""
algorithms: Dict[str, Callable] = {
'Nearest Neighbor': lambda d: nearest_neighbor(d),
'NN + 2-Opt': lambda d: greedy_plus_two_opt(d),
'Simulated Annealing': lambda d: simulated_annealing_tsp(d),
}
# Add exact algorithms only for small instances
exact_algorithms: Dict[str, Callable] = {
'Held-Karp': lambda d: held_karp(d),
'Branch & Bound': lambda d: branch_and_bound_tsp(d),
}
for n in sizes:
print(f"\n{'='*50}")
print(f"Problem size: {n} cities")
print(f"{'='*50}")
distances = generate_random_instance(n)
# Run exact algorithms only for small instances
if n <= 20:
for name, algo in exact_algorithms.items():
start = time.time()
cost, path = algo(distances)
elapsed = time.time() - start
print(f"{name:20} | Cost: {cost:10.2f} | Time: {elapsed:8.4f}s")
for name, algo in algorithms.items():
start = time.time()
cost, path = algo(distances)
elapsed = time.time() - start
print(f"{name:20} | Cost: {cost:10.2f} | Time: {elapsed:8.4f}s")
if __name__ == "__main__":
benchmark_algorithms()
Decision framework:
- Under 20 cities: Use Held-Karp for guaranteed optimal solutions
- 20-50 cities: Branch and bound may work; fall back to SA if it times out
- 50-500 cities: Nearest neighbor + 2-opt offers the best speed/quality tradeoff
- 500+ cities: Simulated annealing or genetic algorithms, with 2-opt for final polish
The greedy + 2-opt combination deserves special mention. It runs in practical O(n²) time and consistently achieves near-optimal results. For most real-world applications—delivery routing, warehouse picking, PCB drilling—this is the right choice. Save the metaheuristics for when you need that extra 2-3% improvement and have the compute budget to find it.