Bitmask DP

Fifteen places in Almaty, from Republic Square and the Green Bazaar to Kok-Tobe and MEGA, and you need to visit them all by road and come back to where you started, driving as little as possible. There are 14!/2 = 43,589,145,600 different tours. Trying them all in Go would take more than ten minutes, judging by the measurement below, while a table finds the best tour, 42.41 kilometres, in hundredths of a second in your browser. Click the points: add up to eighteen, remove some, compare with quick rules.

The rule “always drive to the nearest place not yet visited” covers 48.37 kilometres, 14 % more. It starts well and finishes badly: the last places are scattered across the city, and you have to drive back through half of it to reach them. The dashed “+ 2-opt” is the same rule followed by repairing the route: a piece of the tour is reversed as long as that shortens it. On these fifteen points the repair finds the same best tour, but add a sixteenth, Dostyk Plaza, and it gets stuck at 43.73 kilometres against 42.56. You can only know the exact answer when you have something to compare with.

The idea

What do you need to know about an unfinished route to finish it in the best way? Not the order in which you visited the places, only which places are behind you and where you are standing now. The rest of the way depends on exactly that: which places are left to visit and where you are driving to them from.

The state is a pair “set of visited places, place where we stand”, and the set is written as a number whose bit i is one if place i has been visited. For each state we keep the length of the shortest path from the start that has visited exactly these places and ended at this one. A transition is a trip to a place not yet visited: switch its bit on. There are 2ⁿ sets, not n!: all paths that have visited the same places and stand at the same one merge into one state, and only the shortest of them is remembered.

This is bitmask DP: the same as any dynamic programming, only the state contains a subset. How to work with sets as numbers, uniting them, testing a bit, going through submasks, is told in the article on bit tricks. What matters here is that a mask is an ordinary integer and will do as the index of a table cell.

The Held–Karp table

Number the places from 0 to n − 1 and start at place 0. best[S][j] is the length of the shortest path that left 0, has been to exactly the places of the set S, and stands at j ∈ S. Place 0 is not included in the mask: it is the start.

The answer is the smallest best[all places][j] + d(j, 0): from the last place, return to the start. The tour itself is restored by back pointers: for each state remember which place j it was reached from, and walk from the end, switching bits off.

Each transition switches one bit on, and the mask grows. So the table can be filled in increasing order of masks, or in layers: first all sets of one place, then of two, of three… The layers are shown under the map: the height of a bar is how many states the layer has. Most are in the middle, where the number of sets of s elements, “14 choose 7”, is largest.

The key place in the proof is “the last step does not depend on the piece”. The future of a path depends only on where we are and which places have been taken. If, say, the cost of the next trip depended on the time of day, and the time on the length already driven, the pair “set, place” would not be enough, and the time would have to be added to the state.

Code

Distances between places are by road, by Dijkstra's algorithm from each place; places are attached to the nearest crossroads. Then the table over masks with back pointers, the nearest-neighbour rule, the 2-opt repair, and a check by brute force on the first nine places.

The table for 15 places has 229,376 cells, but half as many real states, 114,688: the place where we stand must be in the mask. There are 745,472 transitions. Compare with the 43,589,145,600 tours of brute force: a difference of tens of thousands of times, and it grows with every place, since brute force multiplies by n while the table only slightly more than doubles. On nine places brute force and the table give the same 27.70 kilometres, as they should.

In Python it is convenient to write out in advance where a mask can step to: the list of free places together with the table rows for the new masks. Then the inner loop has neither a bit test nor an index computation, and the listing runs in a fraction of a second.

What it costs

There are (n − 1) · 2^{n − 2} states, with up to n transitions from each: time O(2ⁿ · n²), memory O(2ⁿ · n). This is an exponential, but incomparably smaller than the n! of brute force.

The time axis is logarithmic. Brute force grows ever steeper: each new point multiplies its time by n, and at 11 points Python takes more than a second and a half. The table grows in a straight line, about five times every two points, and at 16 points takes a tenth of a second in Python. But a straight line on a logarithmic axis is an exponential too. A table for 25 points is 2²⁴ · 24 ≈ 400 million cells, gigabytes of memory. Bitmask DP is a weapon for n up to twenty-odd, and beyond that the travelling salesman problem is solved differently: by branch and bound, linear programming cuts, heuristics with guarantees. That is the subject of a separate article on the travelling salesman problem.

Nobody has yet found an exact algorithm for the general travelling salesman problem, with arbitrary distances, faster than 2ⁿ: papers from 2018 still call the existence of a (2 − ε)ⁿ algorithm an open question. The dynamic programming invented in 1962 is more than sixty years old, and in the worst case it is still the best known.

Sometimes a state needs not one mask but a pass over all its submasks: “split the set into groups”, “cover with sets”. Then there are 3ⁿ transitions; why that many is proved in the article on bit tricks. For twenty elements that is three and a half billion steps, and the estimate is worth making in advance.

A mask as a border: dominoes

A mask does not have to be a set of visited places. In how many ways can an 8 × 8 board be tiled with 1 × 2 dominoes? Lay them cell by cell, row by row. At any moment the tiled part of the board is separated from the empty part by a broken border, and for going on only one thing matters: which cells on this border are already covered by dominoes sticking out of the tiled part. There is one such cell on the border per column, so this is a mask of 8 bits. Press “The profile step by step”: the orange broken line is the border, called the profile.

The transition at each cell is simple. If the cell is already covered from above, switch its bit off and go on. If it is free, lay a domino downwards (the bit stays on, and now it means the cell below us) or to the right (the bit of the right neighbour switches on). A cut-out cell behaves like one already covered, only nothing can cover it. There are 2⁸ = 256 masks, 64 cells, at most two transitions per cell: a few tens of thousands of steps for everything.

For 2 × n strips you get the Fibonacci numbers: the last domino either stands upright, or two lie one above the other horizontally. For 3 × n with odd n the answer is zero, as the number of cells is odd. The 8 × 8 board has 12,988,816 tilings, and this number has been known since 1961: Pieter Kasteleyn and, independently, Harold Temperley with Michael Fisher derived a formula for the number of tilings of a rectangle that involves no search at all. But the formula works only for rectangles, while the profile counts any board, with cut-outs, with holes, with walls.

Cut off two opposite corners. There are no tilings, and the table dutifully finds zero, although the reason is visible without it: diagonally opposite corners are the same colour, leaving 32 cells of one colour and 30 of the other, while each domino covers one of each. This “mutilated chessboard” problem was proposed by the philosopher Max Black in 1946 and made popular by Solomon Golomb and Martin Gardner. But if you cut out two cells of different colours, for example two corners of the same side, tilings always exist (this is Ralph Gomory's theorem): the listing counts 2,436,304 of them.

Where it is used

Join order. A database query joining several tables can be executed in different orders, and the costs of different orders can differ many times over. The PostgreSQL optimiser chooses the order by dynamic programming over sets of tables, in layers, as in the scene: the developer documentation says that the first pass builds joins of two tables, the second of three, the next of four, and so on up to all of them at once. The approach itself goes back to IBM's System R optimiser. The price is known: from twelve tables on, PostgreSQL switches by default to a genetic algorithm (the parameter geqo_threshold is 12), because, as the documentation puts it, with many tables an exhaustive search takes more time than it saves.

Routes with few stops. Visiting a dozen addresses, going round several warehouses, the order of points in a day's plan: where there are few points, the table gives an exact answer in a fraction of a second, and there is no need to guess how far off a heuristic was.

Assignment and packing. Which of n workers gets which job, and how to split things into trips: both are problems below. Assignment also has the polynomial Hungarian algorithm, but masks are simpler and easily take on extra conditions.

Competitions. A limit of “n ≤ 20” in a statement is almost always a hint: 2²⁰ is about a million, and the solution is bitmask DP or meet in the middle.

History

The table for the travelling salesman problem was published in 1962 independently by Richard Bellman, “Dynamic programming treatment of the travelling salesman problem” in the Journal of the ACM, and by Michael Held and Richard Karp, “A dynamic programming approach to sequencing problems” in the SIAM journal. It is called the Held–Karp algorithm, or the Bellman–Held–Karp algorithm. Domino tilings have their own history, older than programs: physicists studied them as the dimer problem, where a domino models a two-atom molecule occupying two neighbouring sites of a lattice, and it was for this problem that Kasteleyn, Temperley and Fisher found their formula in 1961.

Problems

The first is to give sixteen orders to sixteen couriers with the smallest total distance along the streets of Almaty: a mask of orders, and it tells you the courier's number by itself. In the second, move things up in a freight lift in the fewest trips, and one mask is not enough: it needs the weight of the current trip.