Interval DP

The hundred most frequent words of War and Peace in the Russian original, «и», «в», «не», «что», «он» (“and”, “in”, “not”, “that”, “he”) and 95 more, make up 39.7% of all the words of the novel. Put them into a search tree where the words stand in alphabetical order, and search for each as many times as it occurs in the text. How many comparisons will a lookup take on average? If you halve the alphabet, 5.55. If you put the most frequent word at the root, 4.81. And the best tree gives 4.40, and its root is not «и» but «не». The scene builds all three trees, and under them the table that found the best one.

Why not «и», if it occurs more often than any other, 21,079 times? At the root it splits the dictionary very unevenly. To the left of «и» alphabetically there are 32 words, «а», «в», «был», «его»…, and together they occur 46,887 times. To the right there are 67 words, among them «не», «он», «что», «с», «как», «на», and 104,976 uses: more than twice as many. The whole right subtree would go a level deeper, and each of its lookups would get dearer. «Не» splits more evenly, 94,391 uses on the left, «и» itself included, and 70,035 on the right, while «и» ends up on the second level, only one comparison further. The rule “the most frequent at the root” does not think about this and loses 0.41 comparisons on every lookup.

Under the tree is a pyramid. Its bottom row is the intervals of the alphabet of one word, the next of two neighbours, and so on up to the top, where the whole dictionary is. Each block is one cell of the DP table: the best tree for that interval. The colour shows where in the interval its best root stands. Long slanting stripes of one colour are one and the same frequent word that stays the best root for many intervals in a row.

The idea

In a problem about a piece of a sequence, the last (or the first) step divides the piece into two smaller parts somewhere inside: the root of a tree divides an interval of the dictionary, the last multiplication divides a chain of matrices, the first cut divides a log. The parts are solved independently and in the same way. So the state is an interval [i, j], and the transition is trying the place of the cut k inside. Intervals are computed in increasing order of length: everything an interval relies on is shorter than it and already done.

This is interval DP. It differs from ordinary DP over prefixes in that a subproblem is not the beginning of the sequence but any piece of it, and the piece can be cut anywhere. The price for this is an extra factor of n: there are n(n + 1)/2 intervals, and each has up to n ways to be cut.

The best search tree

Words a₁ < a₂ < … < a_n with frequencies w₁, …, w_n. A lookup of a word in a search tree costs as many comparisons as the depth it lies at, counting the root as the first level. The cost of a tree is the sum of w_i · depth(a_i). Let C[i][j] be the cost of the best tree of the words of an interval, and W(i, j) their total frequency. Then

The term W(i, j) appears because the root adds one comparison to the lookup of each word of the interval: all the words of the subtrees go down a level. The sums W are computed in O(1) with prefix sums.

Notice what exactly made the problem an interval problem: a subtree always occupies a solid piece of the alphabet. If the order did not matter, as in the Huffman code, where symbols can be rearranged freely, greedy merging of the two lightest would be enough. Alphabetical order takes this freedom away, and instead of greed the roots have to be tried.

Donald Knuth noticed more in a 1971 paper. The best root of the interval [i, j] lies between the best roots of the intervals [i, j − 1] and [i + 1, j]: if a word is added to an interval on the right, the root can only move to the right. So the root is searched for not over the whole interval but between two already known ones, and the whole table is computed in O(n²) instead of O(n³). When this technique is legitimate and how it generalises is in the article on DP optimisations.

The order of multiplying matrices

Matrix multiplication is associative: (AB)C = A(BC), and the result does not depend on the brackets. But the time does, and a lot. A p × q matrix times a q × r matrix takes p · q · r multiplications of numbers, and in a chain of many matrices different orders give different intermediate sizes. The documentation of the numpy.linalg.multi_dot function has an example of four matrices: 10,000 × 100, 100 × 1000, 1000 × 5 and 5 × 333. Multiplying them from left to right costs more than a billion multiplications of numbers, and the best order 22 million, 48 times fewer.

The problem is about intervals again. The last multiplication divides the chain A_i … A_j into two parts, A_i … A_k and A_k+1 … A_j, and costs p_i−1 · p_k · p_j, where p are the sizes along the chain:

This is exactly what numpy does: in the source of multi_dot the function _multi_dot_matrix_chain_order fills tables of costs and cuts in increasing order of interval length, and a comment says that the implementation “closely follows” the chapter of the textbook by Cormen, Leiserson, Rivest and Stein. The number of orders it chooses from is a Catalan number: 5 for four matrices, 42 for six, more than a billion for twenty. The table goes over all of them without listing a single one: each order is a choice of the last multiplication plus the orders of the two parts.

Switch to “row × … × column”: a row of 1000 numbers, two 1000 × 1000 matrices and a column. If the two big matrices are multiplied first, that is a billion multiplications; if the row is carried through the chain, a couple of million. Hence a simple rule: if the chain has a vector, start multiplying from it, not with the matrices against each other.

Code

The dictionary from the start of the article. The listing counts the words of the novel itself in the sandbox, takes the hundred most frequent, builds three trees and prints the first levels of the best one.

The table is 5050 intervals and 171,700 steps: one step for each “interval, root” pair. The root is «не», under it «и» and «с», further down «в», «как», «он», «что».

The multiplication order: the numpy example and a chain of six matrices from the textbook of Cormen et al., where the best order costs 15,125 multiplications. Together with the best order the listing also finds the worst, with the same table where the minimum is replaced by a maximum.

What it costs

There are n(n + 1)/2 intervals, an interval of length ℓ has ℓ ways to be cut; about n³/6 steps in all, O(n³) time and O(n²) memory for the table. For a hundred words that is 171,700 steps, for a thousand already 167 million.

The slope of trying all roots is 3: twice the words, eight times as long. With Knuth’s window the slope is 2. On four hundred words in Python that is 3.1 seconds against 0.11, almost thirty times faster, and the gap grows. Without such an observation interval DP stays cubic, and n around a thousand is its limit; problems with n up to a hundred or two it solves at once.

Where it works

Linear algebra. numpy.linalg.multi_dot chooses the order of multiplying a chain of arrays with this DP. For three matrices numpy compares the two orders directly, for more it builds the table. If the chain is long, there is a faster algorithm too: T. C. Hu and M. T. Shing found an O(n log n) solution in the early 1980s by reducing the problem to cutting a polygon into triangles.

Biology. An RNA molecule folds onto itself: bases at different ends of the chain form pairs, and the shape of the molecule determines its function. Pairs do not cross: inside a “hairpin” formed by the pair (i, j) all the other pairs lie either wholly inside the interval or outside it. These are intervals again: Ruth Nussinov proposed in the late 1970s a DP that finds the folding with the most pairs in O(n³), and more accurate models count not pairs but energy. The programs of the ViennaRNA package build tables over all intervals of the sequence by the scheme of Michael Zuker and Patrick Stiegler of 1981.

Parsing sentences. The Cocke–Younger–Kasami algorithm checks whether a sentence can be derived from a grammar with the same pyramid: the cell for an interval of the sentence holds which parts of the grammar this interval can be, and it is assembled from two smaller cells.

Problems with cuts and merges. Saw a log, glue neighbouring pieces, put brackets into an expression, cut a polygon into triangles: anything where the last step divides an interval in two and the order of the elements cannot be changed.

History

Best search trees were computed as early as the late 1950s: Edgar Gilbert and Edward Moore solved a close problem about codes that preserve alphabetical order with a cubic DP. Donald Knuth, in the paper “Optimum binary search trees” in the first volume of the journal Acta Informatica (1971), built the table for search trees with frequencies not only of the words but also of the gaps between them, lookups that find nothing, and showed that the roots are monotone and the table can be computed in quadratic time. The matrix chain problem is a classic example from the textbook of Cormen, Leiserson, Rivest and Stein, after whose chapter multi_dot is written.

Problems

The first is sawing a log, where each cut costs the length of the piece: the most direct interval problem. In the second, brackets must be placed so that the expression becomes as large as possible, and one maximum per interval will not be enough.