DP Dynamic programming Technique

Knuth and divide-and-conquer optimisations

Also known as: Knuth optimization · divide and conquer optimization · quadrangle inequality · monotone optimal split · Jenks natural breaks · Monge condition

When the best split in a DP does not move left as the boundary moves right, it need not be searched for anew: divide and conquer computes a layer of a partition into k parts in O(n log n) instead of O(n²), and Knuth's window brings an interval DP down to O(n²). Both rest on the quadrangle inequality for the cost of a piece.

Expert 25 minutes dynamic programming divide and conquer optimization geodata contests

On a choropleth map each town is a dot, and the colour is its population class. Where should the class boundaries go? If the range is divided evenly, almost all 45 thousand towns end up in one class: most have tens of thousands of inhabitants, the largest tens of millions. George Jenks proposed “natural breaks” in 1977: boundaries at which the sum of squared deviations within the classes is smallest. This is a DP: for five classes, trying the boundaries costs five billion evaluations. But the best boundary is monotone, and the same DP fits into three and a half million, which the browser computes in tens of milliseconds.

Five Jenks classes explain 92 % of the spread of population: 44,095 towns up to 527 thousand people, 944 up to 2.2 million, 132 up to 6 million, 44 up to 13 million, and the 12 largest. Quantiles divide the towns evenly by count, and cities of a million drown among towns of a hundred thousand; equal intervals give almost everything to one class.

The idea

Sort the values x1 ≤ … ≤ xn. A class is a run of consecutive values, and its cost w(j, i) is the sum of squared deviations from the mean over the values j + 1 … i; with prefix sums of x and x² it is computed in O(1). The best partition of the first i values into g classes is

bestg[i] = minj < i ( bestg−1[j] + w(j, i) ).

Trying j for all i is n²/2 per layer. But denote by opt(i) the best last split. If the boundary i moves right, the last class should not start further left: opt(i) ≤ opt(i + 1).

If the best split is monotone, compute it first for the middle i = n/2 by trying all j. Then for the left half the splits are no further right than opt(n/2), and for the right half no further left. In each half take the middle again, and the windows narrow. On each level of recursion the windows together cover all j once or twice, there are log n levels, and a layer costs O(n log n).

When the split is monotone

Monotonicity is not a general law. It comes from the quadrangle inequality: w(a, c) + w(b, d) ≤ w(a, d) + w(b, c) for a ≤ b ≤ c ≤ d.

Let w satisfy the quadrangle inequality, and let opt(i) be the largest of the best splits for i. Then opt(i) ≤ opt(i + 1).

Let j′ = opt(i), j = opt(i + 1), and suppose j < j′. The inequality with a = j, b = j′, c = i, d = i + 1 gives w(j, i) + w(j′, i + 1) ≤ w(j, i + 1) + w(j′, i). Since j′ is best for i, f(j′) + w(j′, i) ≤ f(j) + w(j, i), where f = bestg−1. Adding: f(j′) + w(j′, i + 1) ≤ f(j) + w(j, i + 1). So j′ is no worse than j for i + 1 as well, and it is larger, which contradicts the choice of the largest best split.

For sorted values, the cost “sum of squared deviations from the mean over an interval” satisfies the quadrangle inequality.

The sum of squared deviations of a set is the minimum over c of the sum of (x − c)², and the minimum is reached at the mean. Let μ be the mean of the interval [a, d), ν the mean of [b, c) (if it is empty the inequality is obvious: cutting an interval costs no more than not cutting it). Divide [a, d) into three pieces: A = [a, b), B = [b, c), C = [c, d). The right-hand side of the inequality has the terms (x − μ)² over A, B, C and (x − ν)² over B. If μ ≤ ν, give the interval [a, c) the centre μ and [b, d) the centre ν: the left-hand side gets the same terms over A and B, and over C (x − ν)² instead of (x − μ)². The points of C are no smaller than all points of B, hence no smaller than ν ≥ μ, and they are no further from ν than from μ. If ν < μ, give [a, c) the centre ν and [b, d) μ; now the terms over A change, and the points of A are no larger than ν < μ. In both cases the left-hand side with these centres is no larger than the right, and with the best centres all the more so.

Knuth's window

In an interval DP a split is sought for each interval [i, j]: this is how the optimal search tree or the order of merging neighbouring pieces is built. Donald Knuth noticed for the search tree in 1971, and Frances Yao proved in general in 1980: if the cost w satisfies the quadrangle inequality and does not decrease when an interval is widened, the best split of an interval lies between the best splits of the two intervals one shorter: opt(i, j − 1) ≤ opt(i, j) ≤ opt(i + 1, j). Yao's proof is longer, since the same inequality is first carried over to the table of answers itself, and we do not repeat it here.

If for each interval the split is tried only between opt(i, j − 1) and opt(i + 1, j), the whole table is computed in O(n²).

For intervals of one length L the windows cost opt(i + 1, i + L) − opt(i, i + L − 1) + 1. When summed over i, the right end of each term cancels with the left end of the next, and what remains is at most n + (the number of intervals of length L) ≤ 2n. There are n lengths, O(n²) in all.

Code

A layer is the function layer: a stack of intervals [lo, hi] with windows [olo, ohi], without recursion. For five classes on all 45,227 towns Python takes a few seconds. The scene above does the same in the browser for any number of classes.

What it costs

ProblemBrute forceWith monotonicity
Partition into k partsO(kn²)O(kn log n), divide and conquer
Interval DPO(n³)O(n²), Knuth's window
The cost is a minimum of linesO(n²)O(n), the convex hull trick

On a thousand towns brute force in Python is twelve times slower, in Go at four thousand more than two hundred times; on all 45 thousand towns divide and conquer in Go fits into 66 ms.

Where it is easy to go wrong

  • Monotonicity without checking. If the cost does not satisfy the quadrangle inequality, divide and conquer silently gives an answer that is not the best. The sum of absolute deviations from the mean, for example, violates it, already on nine numbers from 0 to 19, while the sum of absolute deviations from the median obeys it. Check on small data against brute force.
  • Unsorted values. The inequality for the sum of squares is proved for sorted values; for an arbitrary order the classes will not be intervals.
  • A window without j < i. The split must be to the left of i: the right edge of the window must be cut down to i − 1.
  • Loss of precision. The sum of squared populations reaches 10¹⁶, and the difference of two close large numbers loses its low digits. This is usually enough for comparing splits, but variants equal up to rounding may be chosen differently in different languages.
  • Recursion in divide and conquer, of depth log n, is safe, and in Knuth's window there is none at all: it is two nested loops.

Where it is used

Choropleth maps: natural breaks are the default classification, for example in ArcGIS. One-dimensional clustering is the same problem as k-means, but on a line it can be solved exactly. Dividing a sequence into k even parts when unevenness is penalised by a square. Interval DPs, such as optimal search trees, merging neighbouring pieces and cutting a stick, with Knuth's window.

History

Donald Knuth in 1971 sped up building the optimal search tree from n³ to n² by noticing the monotonicity of the roots. Frances Yao in 1980 derived it from the quadrangle inequality, a condition named after Gaspard Monge, who noticed it in 1781 in the problem of moving earth. The divide-and-conquer technique for partitions and the stronger SMAWK algorithm (Aggarwal, Klawe, Moran, Shor, Wilber, 1987) use the same monotonicity. George Jenks described natural breaks for choropleth maps in 1977.

Problems

Read “Eugene Onegin” aloud in k evenings as evenly as possible: divide and conquer. And merge neighbouring pieces into one at the least cost: Knuth's window.

Reading evenings Merge neighbours

Data: Places: GeoNames, CC BY 4.0