Dynamic programming

From the Almaty-2 railway station to MEGA Alma-Ata it is 8.7 kilometres in a straight line. How many routes through the streets are there on which each next crossing is closer to MEGA than the previous one? No loops and no steps back. The answer is 1,515,185,802,420,673. Going through such routes at a billion a second would take seventeen and a half days, while the scene below finds this number in milliseconds and does not list a single route. It fills a table: for each crossing, how many routes lead from it to the goal, starting with those nearest the goal. Put the start and finish anywhere and switch what “closer” means.

While the count runs, dots light up in circles around the finish: a dot's colour is how many digits the number of routes from it has. When the table has reached the start, a “river” remains on the map: the thickness of a street shows what share of all routes passes along it. Press “Record of the map”: between two crossings at opposite edges of the map, 14.6 kilometres apart, there are 484,908,221,365,218,608,639,281,781,018 routes. Thirty digits. Enumeration would take more than a thousand ages of the Universe.

Now put the start on Republic Square and the finish by the Central Stadium. There are zero routes: somewhere on the way one has to move away from the stadium, otherwise a block cannot be got round. The table did not make a mistake: it honestly counted that there are no such routes. Switch the rule to “closer by road”, and there are 36 routes: under this rule the shortest path always qualifies.

The idea

A route from a crossing begins with a step to one of the neighbours closer to the goal, and from there it is already a route from that neighbour. So the number of routes from a crossing is the sum of the numbers of its closer neighbours. The same neighbours are needed by many crossings, and counting each anew makes no sense: each number is found once, put in a table, and the table is filled from the goal outwards, so that by the time of counting all the terms are ready.

This is dynamic programming. It has three parts. The state is what one needs to know to solve the rest of the problem: here the crossing you stand on is enough, and how you got there does not matter. The transition is how the answer is assembled from smaller answers: a sum over the neighbours closer to the goal. The order is the order in which to fill the table: by increasing distance to the goal.

States and transitions together form the subproblem graph. On the map it is visible literally: crossings are vertices, and each street turns into an arrow from the far end to the near one. Since the distance to the goal strictly decreases along an arrow, one cannot return along the arrows to where one has been: the graph has no cycles. That is exactly what allows the table to be filled by increasing distance, in topological order. The “A crossing” mode in the scene shows the transition at any place: the green neighbours are closer to the finish, and the crossing's number is the sum of theirs.

The call tree collapses into a table

The same thought can be written as a recursion, and then it is immediately clear why it cannot do without a table. Take a smaller problem: in how many ways can a machine pay out N tenge in coins of 1, 2 and 5 tenge, dropping them one at a time? Order matters here: 2 + 5 and 5 + 2 are different payouts. If the first coin out was c, N − c remains to be paid, so

and terms with a negative argument are simply skipped. f(0) = 1: everything has been paid, and that is exactly one way, paying nothing more. Below are all the calls of this recursion. The width of a rectangle equals its answer, so identical subproblems look the same and are built the same inside.

Without memory the tree for 7 tenge is 62 calls, although there are only eight different questions, from f(0) to f(7). The bottom row shows how many times each was computed: f(1) 15 times, f(0) 26. These are not random numbers. The call f(k) appears in the tree each time the first coins add up to N − k, and that can be done in f(N − k) ways. So the repeats grow as fast as the answers themselves, and the answers grow as a power: each next tenge multiplies them by about 1.7. For 25 tenge and all eight coins, from 1 to 200 tenge, the tree consists of 968,752 calls.

Turn on “with memory”: the first computation of f(k) goes into the table, the other calls take the answer from it, and what is left of the tree is a chain of N + 1 computations. How this is done in code, which keys are needed and where the catches are is in the article on memoisation; here something else matters more. Look at the order in which the bottom row is filled: 0, 1, 2, …, from small amounts to large ones, although the recursion started at the top. Recursion with memory found by itself the order in which each subproblem is ready before it is needed. If this order is known in advance, recursion is not needed at all: a loop from zero upwards fills the same table. This is the bottom-up table, the second way of writing the same DP.

Why the table is right

Reread the first paragraph of the proof: everything rests on the key strictly decreasing. If steps between crossings at the same distance were allowed, a route could go round in a circle, and there would be infinitely many routes. Then the table cannot be filled in any order: a vertex waits for a neighbour, and the neighbour for it. This happens in shortest path problems too: in a graph with cycles DP needs a different state. The Bellman–Ford algorithm is a DP in which the number of edges of the path is added to the vertex, and this number grows with each step, breaking the cycles.

For problems about the best answer the proof is the same, only instead of “all routes are a first step plus a route” one says “the best route is a first step plus the best route from the next vertex”. This property is called optimal substructure, and Bellman stated it as the principle of optimality: whatever the first decision, the remaining decisions must be optimal for the position it led to.

How to devise a DP: five questions

Any DP answers five questions. Let us go through them on the same problem, but with all the coins that can be in a wallet, 1, 2, 5, 10, 20, 50, 100 and 200 tenge, and an amount of 1000 tenge.

1. State. What does one need to know to solve the rest? How much is still to be paid, s. Which coins have already gone out does not matter.
2. Transition. How is the answer assembled from smaller ones? Try which coin goes out next: f(s) is the sum of f(s − c) over coins c ≤ s.
3. Base. Where is the answer known at once? f(0) = 1.
4. Order. In what order to fill so that the terms are ready? By increasing s: each transition reduces the amount.
5. Answer. Which cell of the table answers the question? f(1000).

The result is a number of 233 digits: that many different sequences of coins give 1000 tenge. But usually the question is not about sequences but about sets: how many different handfuls of coins make 1000 tenge if order does not matter. Here the state “how much is left” is not enough, because it cannot tell 2 + 5 from 5 + 2. One also has to remember which coins are still allowed: agree to pay coins from large to small, and then the state is the pair “how much is left, the largest coin allowed”. This trap is taken apart in detail in the course chapter on DP. And the third question, how few coins are needed, lives on the same graph as the first, only instead of the sum over transitions the minimum is taken, and the answer is the length of the shortest path from 1000 to zero.

The table shows the main thing: the question sets the state, not the other way round. If two different positions must be told apart to answer correctly, they must be different states. If there is no need to tell them apart, extra detail only swells the table. For the fewest coins, by the way, the table is not necessary: with tenge coins the greedy way “the largest first” always gives the fewest. But not every set of denominations is so lucky: with coins of 1, 3 and 4, greed pays 6 as 4 + 1 + 1, while 3 + 3 is possible. Where greed is right and where it is not is in the article on greedy algorithms.

Code

The scene from the start of the article. The key is the squared straight-line distance to the goal (an integer, so comparisons are exact) or the distance by road, computed by Dijkstra's algorithm. The DP itself is one function, and it does not care which key: sort the crossings and go through them in order.

The last lines are one pass that answers for all starts at once: the table for the goal MEGA knows the number of routes from each of the 7909 crossings. The most are from a crossing 11.4 kilometres away, a 23-digit number. The record of the map in the scene was found the same way: a pass with each crossing as the goal, 7909 passes, about 45 seconds in Python on the author’s computer.

Coins: recursion without memory and with memory, then the bottom-up table and the three questions from the table above.

968,752 calls against 26: with memory each of the 26 questions, from f(0) to f(25), is solved once. But the program counts 1000 tenge with a table, and with good reason: try orders_memo(500). The first call goes down into orders_memo(499), that one into orders_memo(498), each level adds another frame of the generator expression, and Python stops with RecursionError. A bottom-up table needs no stack depth.

The trap: the number does not fit

Look at the last lines of the output of the routes program:

  79,359,907,956,934,457,242,215 (23 digits)
  in int64: 2,014,951,835,966,190,183

A 64-bit integer holds numbers up to 9,223,372,036,854,775,807, that is 19 digits. An addition beyond the limit is not an error in Go: the number silently “wraps around” modulo 2⁶⁴, and instead of the 23-digit answer a 19-digit one comes out, perfectly plausible to look at. Nothing hints that it is wrong. Python does not know this trouble: its integers grow as much as needed. Go has math/big for such answers, and that is how the listing is written. In olympiads it is done differently: the answer is asked modulo some number, usually 10⁹ + 7. Then each cell holds a remainder, the sum of two remainders is less than 2.1 billion and fits even in a signed 32-bit integer, and remainders can be added at any moment: the remainder of a sum equals the remainder of the sum of remainders.

What it costs

The time of a DP is counted by one formula: the number of states times the number of transitions from each. Each state is processed once, and each of its transitions takes one action.

• Routes through the city. There are V states, the crossings, and as many transitions as streets, O(E) in all. Plus sorting by key, O(V log V). On Almaty that is 7909 crossings and 10,423 additions, one per piece of street.
• Coins. There are N + 1 states, and from each k transitions, the denominations: O(N·k). For sets there are k times more states, but only two transitions from each, and again O(N·k).

One correction: adding thirty-digit numbers is not one action but as many as the numbers have machine words. On routes through the city this is unnoticeable, but where numbers are thousands of digits long, the length of the number must be multiplied into the bound. Memory is a cell per state, and often less: if transitions lead only to neighbouring rows of the table, one or two rows need to be kept.

A measurement: paying out an amount N with all eight coins, by recursion without memory, recursion with memory and a table. The time axis is logarithmic, so an exponential looks like a straight line.

Recursion without memory gets almost three times dearer for every two tenge: that same power 1.7 from the call tree. The other two lines lie at the bottom and hardly grow: by forty tenge it is hundredths of a millisecond. The table is one and a half to two times faster than recursion with memory in Python and still more noticeably in Go: it needs neither function calls nor a dictionary.

Where DP does not work

No order. If subproblems depend on each other in a circle, the table cannot be filled. Sometimes it helps to add to the state a counter that grows with each step, as in Bellman–Ford. Sometimes the order can be found along the way: this is how Dijkstra's algorithm works, taking the nearest vertex from the heap each time, whose answer is already final.

Too many states. To visit all the points of a route once each, the set of points already visited has to go into the state, and there are 2ⁿ sets. For fifteen points this is still bearable, bitmask DP, for a hundred it is not. It can be subtler: in the knapsack problem there are as many states as grams in the capacity, and the table grows with the numbers in the input, not only with how many there are.

No optimal substructure. The longest path without repeated vertices in an ordinary graph cannot be found this way: the continuation of the best path from a vertex may run into vertices already taken by its beginning. To account for this, the set of vertices passed would again have to go into the state.

Where it is used

Communication. The routing protocol RIP, described in RFC 1058 in 1988, is based on the Bellman–Ford algorithm: routers exchange distance tables, and each builds its own from the tables of its neighbours, a DP transition spread over a network. The Viterbi algorithm (1967) recovers a transmitted sequence from a signal distorted by noise, going not through all sequences but through the best paths to each state of the encoder; it is used to decode convolutional codes in cellular, satellite and space communication. More in the article on the Viterbi algorithm.

Texts. TeX breaks a paragraph into lines not greedily, line by line, but for the whole paragraph at once: Donald Knuth and Michael Plass described this in 1981 as finding the best path through a graph of possible break points, without backtracking, by dynamic programming. The “Did you mean…” suggestion and comparing versions of a text rest on edit distance and the longest common subsequence.

Biology and pictures. Aligning genes and proteins: Needleman–Wunsch and Smith–Waterman. Narrowing a photograph without distorting what matters: seam carving, a shortest path from top to bottom through the pixels.

Verses, rabbits and a secretary

One of the oldest rules in the spirit of DP was written down by experts in versification. In Sanskrit prosody syllables are short and long, and a long one lasts twice as long. How many different rhythms of length n short beats can be made? The last syllable is either short, and then before it there is a rhythm of length n − 1, or long, and before it a rhythm of length n − 2. So the number of rhythms of length n is the sum of the two previous ones: the same question “what was the last step” with which any transition begins. This is how Virahanka counted around 700; the original text is lost, but the rule survives in a quotation by the commentator Gopala (around 1135). Around 1150 Hemachandra wrote: the sum of the last and the one before last gives the number of the next metre. The result is the numbers 1, 2, 3, 5, 8, 13, 21, …, the very ones that Leonardo of Pisa, Fibonacci, derived in 1202 in the “Book of Calculation” from the problem of breeding rabbits.

The name appeared in the 1950s. Richard Bellman, a mathematician at the RAND Corporation, worked on multistage decision processes, such as managing stocks, production or rockets, where each decision changes the conditions for the next. In 1954 Bellman set out the theory in a talk at a summer conference of mathematicians, and in 1957 published the book “Dynamic Programming”, which stated the principle of optimality. The word “programming” then meant drawing up plans, not working on a computer. Why Bellman chose exactly this name, and what a defence secretary has to do with it, is told in the course chapter.

Problems

The first is the scene from the start of the article as a problem: routes with no step back, the answer modulo a number. The catch is in the order of counting and in comparing distances. In the second the coins in the wallet run out, and how many coins of each denomination have been spent will have to be added to the state, or a way found to do without it.