Bitmask DP: travelling salesman over subsets

Why this matters

A delivery driver must visit 10 cities and return home. The driver wants the shortest total distance. This is the travelling salesman problem (TSP). The simplest method tries every order of the cities and keeps the shortest route. With 10 cities there are 9 factorial routes, which is 362,880 routes. With 20 cities there are 19 factorial routes, about 1.2 times 10 to the 17th. A computer checking a billion routes per second would still need years.

This lesson builds a faster method. It is still slow as the number of cities grows, but it handles about 20 cities in well under a second. The method is called bitmask DP. Dynamic programming (DP) means solving a big problem by solving smaller pieces once and saving each answer. "Bitmask" means a set of cities is stored as the bits of one whole number. Earlier lessons in this module use DP over numbers and strings. Here the saved answers are indexed by a set of cities and an end city.

The next lesson in this module, on assigning tasks to people, reuses the same bit tricks with a different state.

Sets as bits

A computer stores a whole number as bits. Each bit is either 0 or 1. We give each city one bit. Bit 0 stands for city 0, bit 1 for city 1, and so on. If a bit is 1, that city is in the set.

With 4 cities (numbered 0 to 3), these are some sets:

The number that stands for a set is called a mask. In C++, four operations are enough:

• 1 << c makes the bit for city c. For example, 1 << 2 is 4.
• mask & (1 << c) is not zero when city c is in the set.
• mask | (1 << c) adds city c to the set.
• mask ^ (1 << c) removes city c from the set. It is only safe when city c is already in the set.

The mask for all n cities is (1 << n) - 1. For 4 cities that is 15.

The << symbol shifts bits to the left. It does not mean "less than" here. A set of 20 cities fits in a single int, so checking membership costs one operation.

The state: best route for a set and an end city

We fix city 0 as the start and the finish. Every route begins at city 0. The state is:

dp[mask][last] is the cheapest distance for a route that starts at city 0, visits exactly the cities in mask, and ends at city last.

Two rules keep this small:

• City 0 is always in the mask, because every route starts there.
• last is always in the mask, because the route stands on that city.

Why does only the cheapest route matter? Suppose two routes both visit the cities {0, 1, 3} and both end at city 3. The cities still to visit are the same for both routes. The trip from city 3 onward costs the same no matter how the first routes went. So the route with the lower cost is always the better choice, and the other one can be thrown away.

The rule for moving forward is:

• From state (mask, last), go to an unvisited city next.
• The new state is (mask | (1 << next), next).
• Its cost is the old cost plus the distance from last to next.
• Keep the smallest cost over all ways to reach the same new state.

The start state is the set {0} ending at city 0, with cost 0.

The end is the full set with every city visited, at any end city last. The route must then return to city 0, so the answer is the smallest value of dp[full][last] plus the distance from last back to city 0.

Processing order matters. Adding a city always makes the mask number bigger. So if the program goes through masks from smallest to largest, every state is finished before any state that depends on it is used.

Worked example by hand

Here are the distances for 4 cities. The row is the city you leave, the column is the city you go to.

     0   1   2   3
0    0  10  15  20
1   10   0  35  25
2   15  35   0  30
3   20  25  30   0

Work through the table one set size at a time.

Sets with 2 cities (one trip out from city 0):

• {0, 1} ending at 1: the trip 0 to 1 costs 10.
• {0, 2} ending at 2: the trip 0 to 2 costs 15.
• {0, 3} ending at 3: the trip 0 to 3 costs 20.

Sets with 3 cities:

• {0, 1, 2} ending at 2: come from {0, 1} ending at 1, then 1 to 2. Cost 10 + 35 = 45.
• {0, 1, 2} ending at 1: come from {0, 2} ending at 2, then 2 to 1. Cost 15 + 35 = 50.
• {0, 1, 3} ending at 3: come from {0, 1} ending at 1, then 1 to 3. Cost 10 + 25 = 35.
• {0, 1, 3} ending at 1: come from {0, 3} ending at 3, then 3 to 1. Cost 20 + 25 = 45.
• {0, 2, 3} ending at 2: come from {0, 3} ending at 3, then 3 to 2. Cost 20 + 30 = 50.
• {0, 2, 3} ending at 3: come from {0, 2} ending at 2, then 2 to 3. Cost 15 + 30 = 45.

All 4 cities:

• Ending at 1: the best of 50 + 35 = 85 (from {0, 2, 3} ending at 2) and 45 + 25 = 70 (from {0, 2, 3} ending at 3). The best is 70, the route 0, 2, 3, 1.
• Ending at 2: the best of 45 + 35 = 80 (from {0, 1, 3} ending at 1) and 35 + 30 = 65 (from {0, 1, 3} ending at 3). The best is 65, the route 0, 1, 3, 2.
• Ending at 3: the best of 50 + 25 = 75 (from {0, 1, 2} ending at 1) and 45 + 30 = 75 (from {0, 1, 2} ending at 2). The best is 75.

Close each loop by adding the trip back to city 0:

• Route ending at 1: 70 + 10 = 80.
• Route ending at 2: 65 + 15 = 80.
• Route ending at 3: 75 + 20 = 95.

The shortest tour costs 80. Two tours tie at 80: 0, 2, 3, 1, 0 and 0, 1, 3, 2, 0.

By hand, this took 3 sets of 2 cities, 3 sets of 3 cities (6 states), and 1 full set (3 states). A brute-force search would check 3 factorial equals 6 full routes. The gain is small with 4 cities. The gain grows fast as the number of cities grows.

The program

The program below stores two tables. The dp table holds costs. The parent table holds, for each state, the city visited just before the end city. The parent table lets the program print the route, not only its cost.

The table is one flat list, so the state (mask, last) sits at index mask * n + last. Save this as tsp.cpp:

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <vector>
 
struct Tour {
    int cost = 0;
    std::vector<int> cities;  // starts at city 0; the return to 0 is not listed
};
 
Tour solve_tsp(const std::vector<std::vector<int>>& dist) {
    const int n = static_cast<int>(dist.size());
    const int full = (1 << n) - 1;
    constexpr int kInf = std::numeric_limits<int>::max();
 
    // The state (mask, last) lives at index mask * n + last.
    // dp: cheapest way to leave city 0, visit exactly the cities in mask,
    //     and stand at city last.
    // parent: the city visited just before last on that cheapest way.
    const auto index = [n](int mask, int last) {
        return static_cast<std::size_t>(mask) * static_cast<std::size_t>(n) +
               static_cast<std::size_t>(last);
    };
    std::vector<int> dp(static_cast<std::size_t>(full + 1) * n, kInf);
    std::vector<std::int8_t> parent(static_cast<std::size_t>(full + 1) * n, -1);
 
    dp[index(1, 0)] = 0;
    for (int mask = 1; mask <= full; ++mask) {
        if ((mask & 1) == 0) continue;  // every route starts at city 0
        for (int last = 0; last < n; ++last) {
            if ((mask & (1 << last)) == 0) continue;
            const int cost = dp[index(mask, last)];
            if (cost == kInf) continue;  // this state cannot be reached
            for (int next = 1; next < n; ++next) {
                if (mask & (1 << next)) continue;  // already visited
                const int new_mask = mask | (1 << next);
                const int candidate = cost + dist[last][next];
                if (candidate < dp[index(new_mask, next)]) {
                    dp[index(new_mask, next)] = candidate;
                    parent[index(new_mask, next)] = static_cast<std::int8_t>(last);
                }
            }
        }
    }
 
    // Close the loop: add the trip from the last city back to city 0.
    Tour tour;
    tour.cost = kInf;
    int last_city = -1;
    for (int last = 0; last < n; ++last) {
        const int cost = dp[index(full, last)];
        if (cost == kInf) continue;
        const int total = cost + dist[last][0];
        if (total < tour.cost) {
            tour.cost = total;
            last_city = last;
        }
    }
 
    // Walk the parent links backwards from the last city to city 0.
    int mask = full;
    int city = last_city;
    while (city != -1) {
        tour.cities.push_back(city);
        const int prev = parent[index(mask, city)];
        mask ^= 1 << city;
        city = prev;
    }
    std::reverse(tour.cities.begin(), tour.cities.end());
    return tour;
}
 
// Try every order of cities 1..n-1. Only practical for small n.
int brute_force(const std::vector<std::vector<int>>& dist) {
    std::vector<int> order(dist.size() - 1);
    std::iota(order.begin(), order.end(), 1);
    int best = std::numeric_limits<int>::max();
    do {
        int cost = dist[0][order.front()];
        for (std::size_t i = 1; i < order.size(); ++i) {
            cost += dist[order[i - 1]][order[i]];
        }
        cost += dist[order.back()][0];
        best = std::min(best, cost);
    } while (std::next_permutation(order.begin(), order.end()));
    return best;
}
 
std::vector<std::vector<int>> random_matrix(int n, unsigned seed) {
    std::mt19937 rng(seed);
    std::uniform_int_distribution<int> pick(1, 99);
    std::vector<std::vector<int>> d(n, std::vector<int>(n, 0));
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            d[i][j] = d[j][i] = pick(rng);
        }
    }
    return d;
}
 
void print_tour(const Tour& t) {
    std::cout << "cost " << t.cost << ", route";
    for (int c : t.cities) std::cout << ' ' << c;
    std::cout << " 0\n";
}
 
int main() {
    const std::vector<std::vector<int>> example = {
        {0, 10, 15, 20},
        {10, 0, 35, 25},
        {15, 35, 0, 30},
        {20, 25, 30, 0},
    };
    print_tour(solve_tsp(example));
    std::cout << "brute force: " << brute_force(example) << '\n';
 
    for (int n : {4, 6, 8}) {
        const auto d = random_matrix(n, 100 + n);
        const Tour t = solve_tsp(d);
        const int b = brute_force(d);
        std::cout << "n=" << n << " dp=" << t.cost << " brute=" << b
                  << (t.cost == b ? " match" : " MISMATCH") << '\n';
    }
}

Compile and run it:

g++ -std=c++20 -Wall -o tsp tsp.cpp
./tsp

The compile step prints no warnings. The run prints:

cost 80, route 0 2 3 1 0
brute force: 80
n=4 dp=151 brute=151 match
n=6 dp=217 brute=217 match
n=8 dp=186 brute=186 match

The first line matches the hand calculation: cost 80. The route printed is one of the two tied tours. The program picks it because it checks end cities from 0 upwards and only replaces the best result when a strictly smaller total is found.

The last three lines check the program against brute force on random distance tables with 4, 6 and 8 cities. The random tables are symmetric, so the distance from A to B equals the distance from B to A. The method does not need that, but the test data is simpler to read this way.

Reading the main loop

The table is filled in five steps:

• dp[index(1, 0)] = 0 sets the start state: only city 0 visited, standing on city 0, cost 0.
• The outer loop goes through every mask. The line if ((mask & 1) == 0) continue; skips masks without city 0, because no route can be in such a state.
• The middle loop picks last, the city the route stands on. The check (mask & (1 << last)) == 0 skips cities that are not in the mask.
• The inner loop tries each city next as the next stop. Cities already in the mask are skipped. City 0 is skipped too, since it is always in the mask.
• The if (candidate < ...) line keeps the cheaper route. When a cheaper route is found, the program stores its cost and saves last in the parent table.

The index uses std::size_t, the unsigned type C++ uses for sizes and positions in memory. It is the usual type for indexes into a std::vector.

Following the parent links

The route is rebuilt from the end. Start at the best end city, look up its parent, remove that city from the mask, and repeat. For the 4-city example the program keeps end city 1, because the tie at city 2 does not replace it. Its parent is 3, so the route ends 3, 1. Remove city 1 from the mask, which leaves {0, 2, 3}. The parent of city 3 in that set is 2. The parent of city 2 in {0, 2} is 0. The loop stops when the parent is -1, which happens at the start state {0}. The list is reversed at the end, so it reads 0, 2, 3, 1.

The program does not need the parent table to find the cost. Only the route printing uses it. Storing it costs one byte per state, as explained below.

A second example: the nearest city is not always the answer

A common first idea is the nearest-neighbour rule. From the current city, always go to the closest city that has not been visited yet. It is quick, and it often gives a good route. It can also miss the best one. Here are 5 cities:

     0   1   2   3   4
0    0   2   9  10   7
1    2   0   6   4   3
2    9   6   0   8   5
3   10   4   8   0   6
4    7   3   5   6   0

Run the nearest-neighbour rule from city 0. The closest city to 0 is city 1, at distance 2. From city 1 the unvisited cities are 2, 3 and 4, at distances 6, 4 and 3, so it goes to city 4. From city 4 the choices are city 2 (distance 5) and city 3 (distance 6), so it goes to city 2. Then city 3, at distance 8, and back to city 0, at distance 10. The total is 2 + 3 + 5 + 8 + 10 = 28.

The DP program, the brute-force check and a short nearest-neighbour loop were run on the same table:

nearest neighbour cost 28, route 0 1 4 2 3 0
cost 26, route 0 4 2 3 1 0
brute force 26

The best tour costs 26, which is 2 less than the nearest-neighbour tour. The best tour starts with the longer first step from 0 to 4 (distance 7), not the shortest one to city 1 (distance 2). The table keeps the state for the cities {0, 4} standing on city 4, because that state's cost is the cheapest way to reach it. The nearest-neighbour rule throws that state away at the first step. The DP keeps every state until the end and only discards states that are more expensive than another route to the same set and end city.

This table has 2^5 times 5, which is 160 states. The DP checks all of them in a fixed order, so it cannot be misled by an early choice.

How much work the table takes

For n cities the table has 2^n * n states. Each state is visited once, and each visit tries up to n next cities. The total work is at most 2^n * n * n simple steps.

The timing below uses a random table of 20 cities, built with the same random_matrix function and seed 7. It was compiled with -O2, the optimisation flag that makes the compiler produce faster code:

n=20 cost 273 cities 20
seconds 0.344622

That run took about 0.3 seconds in this environment, and the time will differ on your machine. A build without optimisation is slower, but the same order of magnitude is typical. The numbers below are counts for n = 20:

• States: 2^20 times 20 is 20,971,520.
• The dp table stores one int (4 bytes) per state: about 84 MB.
• The parent table stores one byte (std::int8_t) per state: about 21 MB.
• Inner-loop steps: at most 2^20 times 20 times 20, which is about 420 million.

Each extra city doubles the number of sets. At 25 cities the table would have 2^25 times 25, about 839 million states. Each state would need several bytes, which is several gigabytes. This is why the method stops being practical around 20 cities.

Common mistakes

Forgetting the return trip

The loop must end at city 0. If the program takes the minimum of dp[full][last] without adding dist[last][0], it finds the cheapest open path instead of the cheapest tour. For the 4-city example, the cheapest open path ends at city 2 with cost 65, not the tour cost 80. The answer is too small, and the printed route does not return home.

Adding to an unreachable state

Some states cannot be reached. The program marks them with kInf, the largest int value. Adding a distance to that value makes it overflow. Signed overflow is undefined behaviour in C++, so the program may print anything. Compiling with the undefined behaviour sanitizer shows the problem directly. This is a small test program written only for this lesson:

#include <iostream>
#include <limits>
int main() {
    const int inf = std::numeric_limits<int>::max();
    int dp = inf;
    int cost = 10;
    std::cout << dp + cost << "\n";
}

g++ -std=c++20 -Wall -fsanitize=undefined -o inf_add inf_add.cpp
./inf_add

The output is:

inf_add.cpp:7:31: runtime error: signed integer overflow: 2147483647 + 10 cannot be represented in type 'int'
-2147483639

The fix is the line if (cost == kInf) continue; in the program. A state is only extended when its cost is a real number.

Missing brackets around a comparison

In C++, == binds more tightly than &. So mask & (1 << last) == 0 does not test whether the bit is clear. It computes mask & ((1 << last) == 0), which is almost always 0. The compiler warns about it with -Wall:

precedence.cpp:5:28: warning: suggest parentheses around comparison in operand of ‘&’ [-Wparentheses]
    5 |     if (mask & (1 << last) == 0) {

The fix is if ((mask & (1 << last)) == 0), which is the form the program uses.

Treating the table as a 2D array

The table is a flat std::vector, not a grid of rows. Writing dp[mask][last] on a std::vector<int> is a compile error. This three-line file, saved as bad_index.cpp, fails like this:

#include <vector>
int at(std::vector<int>& dp, int mask, int last) {
    return dp[mask][last];
}

The compiler reports:

bad_index.cpp:3:20: error: invalid types ‘__gnu_cxx::__alloc_traits<std::allocator<int>, int>::value_type {aka int}[int]’ for array subscript

The compiler is telling you that dp[mask] is already an int, which cannot be indexed again. Use dp[index(mask, last)] instead, as the program does.

Going above 20 cities

Setting n to 25 does not crash at once. The program asks for about 839 million states. On a machine with a few gigabytes of memory, the std::vector allocation fails with std::bad_alloc or the system stops the program. Check the size before running. For a number of cities around 20 or more, use a method for approximate answers instead.

Practice

Each exercise is followed by its answer.

Exercise 1

For the 4-city example in this lesson, what are dp[{0, 1}][1] and dp[{0, 1, 2}][2]?

Show the answer

dp[{0, 1}][1] is 10, the trip from 0 to 1. dp[{0, 1, 2}][2] is 45, the trip 0, 1, 2, which is 10 + 35.

Exercise 2

How many states does the table have for n = 10? Use the formula 2^n * n.

Show the answer

2^10 is 1,024. Times 10 cities gives 10,240 states.

Exercise 3

Three cities have these distances. Find the shortest tour that starts and ends at city 0.

     0   1   2
0    0   1   2
1    1   0   3
2    2   3   0

Show the answer

There are two tours. Going 0, 1, 2, 0 costs 1 + 3 + 2 = 6. Going 0, 2, 1, 0 costs 2 + 3 + 1 = 6. The answer is 6.

Exercise 4

Why does the loop over masks go from 1 to full in increasing order?

Show the answer

Adding a city to a mask always makes the mask number bigger. So when the loop reaches a mask, every smaller mask that can lead to it has already been finished. Its cost is then final and can be used to extend it.

Exercise 5

How many routes does the brute-force function check for 10 cities? How many does it check for 11 cities?

Show the answer

For 10 cities it checks 9 factorial, which is 362,880 routes. For 11 cities it checks 10 factorial, which is 3,628,800 routes. The number of routes grows by a factor of 10 with each extra city, which is why the table method is needed.

Summary

• A set of cities is stored as a mask, with one bit per city.
• The state is the cheapest route that starts at city 0, visits exactly the cities in the mask, and ends at the last city.
• Only the cheapest route for each (set, end city) pair needs to be kept.
• The answer is the smallest closing cost over all end cities, which includes the trip back to city 0.
• The method stores about 2^n times n states, so it handles about 20 cities in well under a second but not much more.