Dynamic Programming in Practice: Optimal Data Segmentation

I have been reviewing dynamic programming recently, and today I saw a nice application of dynamic programming used in Astrophysics: Dynamic Programming in Python: Bayesian Blocks.

The following is the problem in astrophysics phrased in a purely layman’s way: given a set of data points in the plane, each point is a pair — (time point, observation value), we would like to extract as much information as possible from these observations — which roughly speaking, involves computing a histogram of the observations. The difficulty here is determining the bin size: if the bin size is too large, we loose important information; however, at the same time, we would like to suppress the noise in the observations. It seems that one approach used in the past is by adjusting the bin size manually by a trial-and-error fashion: if the histogram is too jagged, then decrease the bin size; and if it is too smooth, increase the bin size.

The novel idea in this paper is that one first decides a fitness function/cost function over partitions of data using a Bayesian likelihood framework, and then tries to compute an optimum partition of the data points according to this function.

Here is the same problem phrased in the language of computer science:

Let be a set of points in the plane. A block of points in is a set of points consecutive in -axis. A partition of into block is specified by a sequence of blocks , where the blocks are pairwise disjoint, and . Let be the fit of the block . The goal is to find an optimum partition , maximizing (the value is not given, it is also determined by the algorithm).

First renumber the points in by their -coordinates. The recursion formula is straightforward: let denotes the optimum partition of the points in (the leftmost points in ), for and let . Then . There are totally entries to compute, and we spend time for the entry, so the overall time complexity is . Note that we assume that the function can be obtained in time for a block, which usually is true (by possibly preprocessing).

An almost identical problem is segmented least square. Given a set of points in the plane, we would like to partition into blocks, and fit each block with a line. The overall cost the partition is the summation of the costs from each block.

This problem is a generalization of linear regression: in linear regression, the cost of fitting a point set of size to a line is
, where . The “best” coefficients and , that is, minimizing , has closed forms:
,

.

However, sometimes the data points lie on more than one lines. Hence we would like to discover this fact (for example, we might be interested in (1) the number of lines; (2) the points where the line changes direction — these kind of points are called change points in the previous Bayesian block partition problem) and find the optimum cost. A human being can easily determine this sometime, but naturally there are situations where the number of lines we need to “best fit” the point set is quite unclear.

We now formulate this problem into an optimization problem. Given a partition of the point set into segment/block, where each segment consists of a subset of contiguous in -axis, , let be the optimum line fitting , the cost of the partition is
,
where is a positive real number, which encodes the penalty for partitioning into parts (without this term, we could fit each single point with a line, causing a trivial and meaningless optimum solution of cost 0). The goal is to find an optimum partition minimizing .

The dynamic programming approach is the following. First, we renumber the points in according their -coordinate (in ascending order). The recursion formula is the following: define , and define to be an optimum partition of the first points of . The recursion formula is almost identical to the Bayesian Block partition problem:
, for , where is the optimal cost to fit the points in with a single line (we are slightly abusing the notation here).

The time complexity is , following the same reasoning. One thing to note is that we need to do some preprocessing on , so that we can quickly compute in time.

Dynamic Programming Practice (4)

Shuffle of strings A shuffle of two strings and is formed by interspersing the characters into a new string, keeping the characters of and in the same order.

Example: the string “prodgyrnamammiincg” is a shuffle of “dynamic” and “programming”.

Given three strings , and , describe and analyze an algorithm to determine whether is a shuffle of and .

This problem is from these lecture notes.

This is a dynamic programming problem. The recursion formula is straightforward: define to be true if is a shuffle of and , and false otherwise, for and .

(1) . (An empty string is a shuffle of two empty strings).
(2) If , then is true if and only if the string is a prefix of , for .
(3) If , then is true if and only if the string is a prefix of the string , for .
(4) Now assume that both and are non-empty strings, then:

• If , then equals to if , and if ; otherwise it is false.
• If , then .

is the final answer.

The time complexity of the algorithm is , the space complexity is . Here is the result of the execution of the algorithm on the example .

Note that we do not need to explicitly maintain the whole table in the implementation: we fill up the table column by column from left to right, and for each column, the entries are filled from north to south. We only need to maintain one column of the table (the “sliding window” method).

One thing to note in this table is that most of the entries seem to be false, so is there any linear time (that is, ) algorithm for this problem?

A probability puzzle: loaded dice

I came across this question on interview street:

Random number generator
There is an ideal random number generator, which given a positive integer M can generate any real number between 0 to M with equal probability.
Suppose we generate 2 numbers x and y via the generator by giving it 2 positive integers A and B, what’s the probability that x + y is less than C? where C is a positive integer.

Input Format
The first line of the input is an integer N, the number of test cases.
N lines follow. Each line contains 3 positive integers A, B and C.
All the integers are no larger than 10000.

Output Format
For each output, output a fraction that indicates the probability. The greatest common divisor of each pair of numerator and denominator should be 1.

Input
3
1 1 1
1 1 2
1 1 3

Output
1/2
1/1
1/1

One way to consider this problem is by generating functions. Suppose the distribution is not uniform, is there any efficient way to compute the probability that the sum of the two randomly generated numbers is no greater than ?

Suppose for the first random number generator, the probability of that the randomly generated number, , equals to , is , for , and similarly, for the second random generator, the probability that the randomly generated number is is , Define two polynomials and . Then consider the multiplication of :
,
The coefficient in the right hand side of the above equation is precisely the probability of . So the probability that is the summation of the first terms. In the problem of interview street, since the distribution is uniform, each is , and each is . So we can easily obtain the probability that .

The probability puzzle (using generating functions) is the following:
Loaded dice: No matter what two loaded dice we have, it cannot happen that each of the sums 2,3, …, 12 comes up with the same probability.
(I first came across this puzzle from the book “The art of mathematics”).

We can reason like this: let denotes the probability that the result of the first dice is , for . Similarly, let denotes the probability that the result of the second dice is , . Then the probability that the sum of the two results is is (the first dice and second dice must both output 1). The probability that the sum is 12 is (both dices must output 6). Then consider the probability that the sum is 7: it is at least (the first dice outputs 1, and the second dice outputs 6, or the reverse). But .

This reasoning works, but the book provides a much better proof using generating functions. As before, let and let . Now if the sums 2,3,…,12 comes with the same probability, 1/11, then we have
. Now since both and are polynomials of degree 5, has real roots; however, the right hand side is a cyclotomic polynomial which does not have any real root. Thus the identity cannot hold.

Here is another application of generating functions for computing probabilities.

Counting heads Given integers and , along with , you want to determine the probability of obtaining exactly heads when biased coins are tossed independently at random, where is the probability that the coin comes up head. Assume that you can multiply and add two numbers in in time.

(This problem is from this book.)

This problem can be solved using dynamic programming: Let be the indicator random variable for the event that the result of tossing coin is head, for . Notice that
.

The time complexity is . However, if we consider the approach generating functions, we can get a faster algorithm using fast fourier transform:
Let , for . Then consider the product of the polynomials of degree 1: . The probability of obtaining heads is the coefficient . Therefore, the problem is reduced to how fast we can multiply polynomials. Using fast fourier transform, we can multiply two polynomials of degree at most in time . Hence, we can use divide-and-conquer to multiply these polynomials:
.
Let denotes the time to multiply polynomials of degree 1. Then we have the following recursion:
,
as in the “conquer” step, we are multiplying two polynomials of degree . Hence .

Coin Exchange: Dynamic Programming Practice (3)

Coin Exchange: Given a set of denominations , and a target value , find the smallest number of denominations whose sum equals to . For example, if the denominations are , and the target is 122, then the optimum way to pick the denominations is , so the number of bills we need is 5.

This is a very common dynamic programming problem. To simply things a bit, we assume that the smallest denomination is 1, which ensures that for any target value, there is at least one way to get that target (a sequence of 1’s). W.l.o.g., assume that . The recursion formula is the following: define to be the smallest number of bills we need to get , with the constraint that the largest denomination we can use is , for . The base cases are: , for ; for . For and , if , then , where is the largest index satisfying ; otherwise, .

Time Complexity/Space complexity: there are several ways to implement the dynamic programming algorithm. If we use the filling-up table approach to implement the algorithm, the time complexity is usually , where is the number of denominations and is the target value. So it is a psedu-polynomial time complexity algorithm. The space complexity is . In the code below, I used “recursion + memoization” approach (the reason is given here).

#include <vector>
#include <algorithm>
#include <unordered_map>
#include <iostream>
using namespace std;

size_t hasher(const pair<int, int> &p)
{
    return hash<int>()(p.first) + hash<int>()(p.second);
}

class CoinExchange {
public:
    CoinExchange(const vector<int> &Denominations): Denom {Denominations}, cache(Denom.size(), hasher) {sort(Denom.begin(), Denom.end());}
    int operator()(int k);
    int MinCoin(int target, int largest);
private:
    vector<int> Denom;
    unordered_map<pair<int, int>, int, decltype(hasher) *> cache;
};

int CoinExchange::MinCoin(int target, int largest)
{
    if (target <= 0)
        return {};

    if (!largest)
        return target;

    if (target < Denom[largest]) {
        largest = upper_bound(Denom.begin(), Denom.end(), target) - Denom.begin() - 1;
        return MinCoin(target, largest);
    }

    if (cache.find({target, largest}) != cache.end())
        return cache[{target, largest}];
 return cache[{target, largest}] = min(1 + MinCoin(target - Denom[largest], largest), MinCoin(target, largest - 1));
}

int CoinExchange::operator()(int target)
{
    return MinCoin(target, Denom.size()-1);
}

int main()
{
    vector<int> denominations {1,4,7,13,28,52,91,365};
    int target;
    cin >> target;
    CoinExchange coinexchange {denominations};
    cout << coinexchange(target) << endl;
}

Update:
See here for a very well-explained article for solving the coin exchange problem in time.

Dynamic Programming Practice (2)

Maximum set of non-intersecting lines Consider two horizontal lines and in the plane. Suppose we are given the -coordinates of distinct points on and distinct points on . Design an algorithm to compute the largest set of non-intersecting line segments satisfying to the following restrictions:

(a) Each segment in connects some point to the corresponding point .
(b) No two segments in intersect.

This problem is from the chapter of dynamic programming of these lecture notes.

It is quite easy to derive a recursion formula for this problem: let denotes pairs of segments corresponding to , . Then the maximum set of non-intersecting segments of is , where is the set of segments in intersecting . However, in order to use the dynamic programming approach, it is not clear how the subproblems overlap and in what order should we evaluate them.

An alternative approach is to formulate the problem as a graph problem: each segment correspond to a vertex in the graph, and two vertices are connected by an edge if the segment corresponding to them intersect. Then the problem is equivalent to finding a maximum independent set in the graph. However, only when the graph instance has some special properties, we can compute the maximum independent set in polynomial time easily (via dynamic programming); for example, when the underlying graph is a tree. But the graph representing the intersections of segments clearly does not have such a property.

So how to use dynamic programming to solve this problem? The key observation is the following: suppose we renumber the indices of , so that they are ordered from left to right, and we also renumber the corresponding indices of ‘s (so when we change the index of some to , we also change to ). Then we looked at the indices of ‘s. The key observation is that a maximum non-intersecting set of segments correspond to the longest increasing subsequence. The renumber step can be done in , and it takes time to find the longest increasing subsequence. So the overall time complexity is .

Largerst k elements in a stream

Given a sequence of numbers, whose length is unknown, find the largest elements.

This problem is an interesting variant of the largest $k$ elements problem, which is usually phrased as the following: given a sequence of numbers, and a number , find the largest elements. There are two ways to solve this problem: (1) use the selection algorithm to find the element of rank in linear time, and then do another scan to find all elements no larger than the rank element. The time complexity is , space complexity . (2) we can reduce the space complexity by maintaining a min heap of size . For each element after the element, if it is smaller than the top element in the heap, discard it, otherwise, pop out the current smallest element in the heap, and push this new element into the heap. The time complexity is , and the space complexity is .

How to achieve the time complexity while still using space? The solution is the following: use an array of size . If the array is not filled up, add new element into the array. If the array is filled up, use the median selection algorithm to find the median of current element, and discard the elements which are greater than the median. This takes time. The array will be filled up times, so the time complexity is .

I came across this solution from Algorithms Notes.

Breaking a string into words

You are given a string of characters , which you believe to be a corrupted text document in which all punctuation has vanished (so that it looks something like “itwasthebestoftime…”). You wish to reconstruct the document using a dictionary, which is available in the form of an Boolean function : for any string , returns true if is a valid word, and false otherwise.
Given a dynamic programming algorithm that determines whether the string can be reconstituted as a sequence of valid words. The running time should be at most , assuming calls to takes unit time.

The question is an exercise from “Algorithms” by S. Dasgupta, C.H. Papadimitriou, and U.V. Vazirani. It is a typical dynamic programming problem.

The recursion is straightforward: define to be a function over strings, where (the empty string is considered as a valid word in dictionary, and thus a valid sentence vacuously); for ,
.
(The string can only be broken into valid words if and only if there exists some index , such that the string is a word in the dictionary, and the remaining string can be broken into a sequence of valid words). We can slightly simplify the notation by dropping the end index of the string:

for

Then we only need to compute .

The code below is “recursion + memoization” (note that “memoization” should not be confused with “memorization”). In the past I have always thought that “recursion + memoization” is roughly equivalent to iteration (since it is also very straightforward to implement the recursion in an iterative fashion). I have always slightly preferred the iterative fashion over recursion calls, since the latter involve extra bookings such as  pushing parameters on the stack and returning value, etc., though it is indeed much easier to code, and mathematically more concise than the iterative form. But the authors in the book point out an interesting situation where recursion + memoization might be better:

In some cases, though, memoization pays off. Here’s why: dynamic programming automatically solves every subproblem that could conceivably be needed, while memoization only ends up solving the ones that are actually used. …..The memoization recursive algorithm will never look at these extraneous table entries.

bool validsentence(const string &s, size_t i, unordered_map<size_t, bool> &cache, const Dictionary &dictionary)
{
if (i == s.size())
return true;

if (cache.find(i) != cache.end())
return cache[i];

bool res {};
for (size_t j = i+1; j <= s.size(); ++j)
res = res || dictionary.isword(string(s.cbegin() + i, s.cbegin() + j)) && validsentence(s, j, cache, dictionary);
cache[i] = res;
return res;
}

bool validsentence(const string &s, const Dictionary &dict)
{
unordered_map<size_t, bool> cache;
return validsentence(s, 0, cache, dict);
}

Longest subarray whose sum <= k

Longest subarray whose sum
Given an array of numbers and a key , find the longest subarray of for which the subarray sum is less than or equal to .

This problem is from the book “Elements of Programming Interviews”.

This problem has a rather easy quadratic time complexity algorithm: check each pair of indices , where to see if the sum of the subarray is no greater than , and return the pair with the largest difference . To avoid repeatedly computing the subarray sum, one can first preprocess the array to get an array of partial sums, where , for , and by definition, . Then the subarray sum can be obtained as . The preprocessing takes linear time and linear space.

In the following we describe an algorithm with time complexity .

Algorithm LONGESTSUBARRAY(, )
Input: An array and a real number
Output: A pair of indices , maximizing , subject to
1. Compute the array of partial sums, where is the prefix sum of the first numbers in . .
2.
3.
4.
5. for down to 0
    if
        
        
6.
7. for up to
    if
        Use binary search to look for the rightmost element in such that . If such an element exists, update to if
        
8. return

Correctness of the Algorithm
The key is to observe the following two facts.
Claim: (1) If two indices satisfies that , then cannot appear in an optimum solution; (2) If two indices satisfies that , then cannot appear in an optimum solution.
Proof: For any index , , note that the subarray sum of is less than the subarray sum of , and the length of the subarray $A[i..r’]$ is greater than the length of the subarray (although both length might be negative, but the statement still holds). Therefore, is preferable over . This is also the reason why we compute the strictly increasing sequence of .
The second statement follows a similar reasoning.

Time Complexity and Space Complexity
Step 1, computing the prefix sums, takes time. Step 4, computing the array , takes time. Step 7 takes time as each iteration takes time.

Task Scheduling and Balanced Binary Search Tree

The Task Scheduling problem is from Hackerrank.

Task Scheduling
Given a list of tasks , where the task has deadline and duration , for each , , compute the minimum maximum lateness of the first tasks.

First we introduce some notations which will be used in the following discussion. Given a set tasks, a schedule of the tasks is a sequence , where each is in . In other words, the sequence, is a permutation of . For example, a schedule for , could be (we first process task , then , lastly, ). For a schedule , the start time for , , is (the summation of durations of all the tasks before task in the schedule), and the finish time is (the start time for plus the time needed to process ). The lateness of , for in this schedule is the difference between its finish time and its deadline, . Here is an example (the input is from Hackerrank). Suppose we have 5 tasks, . The first task has a deadline at 2, and we need 2 time units to process ; similarly, the second task has a deadline at 1, and its duration is 1, and so on. A valid schedule for these 5 tasks might be, . The start time , for is 0, and the finish time is 1, the lateness is ; the start time for task is 1 (which is the finish time of the task just precedes it), and the finish time is , so it has lateness . Similarly, we can compute the start time/finish time/lateness for the rest of the tasks: , , ; , , ; , , . So the maximum lateness, is 3, which is the lateness of .

A Greedy Algorithm to Process a Set of Tasks
If we are given tasks, and our goal is to find an optimum schedule for the tasks, we can solve this problem using the greedy algorithm: sort the tasks by deadlines in non-decreasing order, and that order is the schedule minimizing the lateness. The time complexity is .

The algorithm, which we will refer to , is the following:

Input: a set of tasks, each task is specified by
Output: a schedule , where , and the maximum lateness .
1. Sort the tasks by the deadlines, and renumbering these tasks. Let be the new sequence of tasks.
2.
3.
4. for each in 

if  ()

5. Output and

The time complexity of the above algorithm is , as the first step, sorting, takes . The rest takes .

An Incremental Algorithm to Process sets of tasks
Let us now consider the Task Scheduling problem from Hackerrank. Here we need to compute the maximum lateness for sets of tasks, that is, the first tasks, for each . For notational convenience, in the following, we use to denote the first tasks. Note that . The possible solutions might be:

• The most naive method. Repeat for each , for . The time complexity is , which is .
• The relatively less naive method. Note that we don’t need to repeatedly do the sorting step: once we have the optimum schedule for , we can use binary search to find the position to insert , which takes . But we still have to update the lateness. How can we do this? Since all the tasks after the newly inserted task is delayed further by time units, their lateness might increase. So if we scan all tasks to find the new maximum lateness, this would take time. So the overall time complexity is , which is , which is .

In the following, we describe an algorithm that takes time and space. Note that the main difficulty is to update the lateness, which in turn involves two subproblems: (1) compute the lateness of the newly inserted task , so that we can compute its latesness , and (2) compute the maximum lateness of all the tasks that are pushed time units further to the right along the time line, because of the insertion of .

Let us first consider the first subproblem, which can be solved in logarithmic time in many ways. The easiest method to think of, is probably to use a balanced binary search tree (but note that we will use a better approach — a simpler one, although with the same time complexity). Given a list of tasks, we can build a vanilla balanced binary search tree of the tasks, where the key is the deadline (this would allow us to find the correct position for the new task in time, as the tree height of a balanced binary search tree with leaves is ). We can augment in the following way: store in each node not only and , but also the summation of the durations of all its descendents. This is the time that is needed to process all tasks in the subtree rooted at this node. So each node in has three values, , and .  When we insert a new task into , we can update this field of all ancestors of and compute the finish time of by the time we reach the bottom of simultaneously: initialize with 0. We will accumulate the time needed to process all tasks before in . Start from the root of , for each node we touched, add to . If the deadline of the new task, is greater than or equal to , descend to the right, and add the of the left child of and to , as the tasks in the left subtree of , and itself, are processed before in the schedule; otherwise, is less than , in this case we descend to the left child, and do not add any value to . When we reached the bottom of , we insert a new node , for , , , and (as has no descendent other than itself). is now the partial sum , where the summation is over all the predecessors of , in other words, all the tasks processed before in the current schedule. With the finish time for , we can compute the lateness of in the current schedule as . The time complexity is linear in the height of , as we only spend time per node as we walk down the tree. Since is balanced, the time complexity is thus .

Let us now consider the second subproblem. How can we compute the maximum lateness of current schedule? Suppose the maximum lateness occurs after , then we can easily get the new maximum lateness: consider any task in . If this task is before in the new schedule, then the insertion of has no effect on its lateness. The insertion of causes a shift for all the tasks in preceded by by time units. Therefore, their maximum lateness, , is now increased by . Hence, the new maximum lateness is . What about the case where the maximum lateness occurs before ? Because the lateness of all the tasks preceded by has increased, it is possible that one of them assumes the maximum lateness. In the worst case, if is inserted to the very front of the schedule, there would be tasks after , even if we do not have to recompute the lateness of each of them,  it is quite inefficient to get the maximum lateness by scanning each of them.

There are two optimizations we could make here. The first one is an observation of the lateness, the other is the use of Fenwick Tree. Let us first consider the following question: in a schedule of tasks, , ordered by their deadlines, which of these tasks might assume the maximum lateness after an insertion of a new task? Suppose we have a sequence of lateness , where is the lateness of task in the schedule. We make the following claim:

Claim: Consider a pair of tasks and , where (that is, task is precedes task in the schedule). If , then task can never have a larger lateness that task .

Proof: The reason is that (1) if the new task is inserted after , both and are not affected; their lateness is the same after the insertion; (2) if the new task is inserted after and before , the lateness of , remains the same, while increased; (3) if the new task is inserted before , then both and is shifted further to the right, by an amount of — the duration of the new task, hence if is no smaller than , it will remain so after the insertion.

This observation allows us to store a sub-sequence of the lateness sequence , which is , where and is a sub-sequence of . Let us denote this list of lateness as . Now we can compute the maximum lateness after the insertion of the new task: find the position where the lateness of would occur in . Suppose and are the two indices. The maximum lateness is the maximum of the lateness before (as all tasks before and including precedes task , so their lateness remain the same after the insertion), the lateness of , and the maximum lateness of all the tasks after and including plus , which is the duration of the new task. Since the elements in is in strictly decreasing order, we can easily the get the maximums of the two parts in : the very first element of each part. Therefore, determining the maximum lateness is . However, we need to update the list : namely, we need to add a value of to each entry after and including , and possibly inserting (the lateness of the new task) before , and remove all the elements in the left part of in order to maintain the strictly increasing order in . (By the way, even if we don’t optimize further, and just naively update , the program already can pass all the tests on Hackerrank).

In the worst case, the size of still could be in the order of the number of all the tasks. So the overall time complexity to update could be large (for example, if is inserted before all the tasks in and each task in the optimum schedule of happens to have a decreasing lateness). While practically fast enough, the time complexity is still not yet. To avoid many updates in after an insertion, we use a simple trick: instead of explicitly storing the lateness in the list, we store the deadline. That is, if appears in , in the same position, we would store , which is the deadline of task . The deadline would serve as an id for each task. Using the deadline, we can find the finish time for the task with deadline and hence compute its lateness in the current schedule in time .

Now in an insertion, there might be several queries on the balanced binary search tree for the finish times, namely, in the set of  tasks precedes in the schedule, which also appears in , we have to find those tasks whose lateness are smaller than the current maximum lateness and throw them out from . This can be done by scanning the elements in reversely from the point where is inserted, till we find a task that precedes with a larger lateness, or the list is exhausted. Since each task can only be thrown away once (if it is evicted from , it would never appear in ), therefore, the overall costs of all insertions is .

After Thoughts/Alternatives
Instead of using balanced binary search tree, we could also use Fenwick Tree. Fenwick tree allows one to update an array and query partial sums on the array efficiently: each operation takes time. In some sense, it is equivalent to balanced binary search tree. But it is much easier to implement . Although conceptually not too hard, it still takes much more lines to code an augmented balanced binary search tree than Fenwick Tree. There is a balanced binary search tree in STL, but currently I don’t know any good to augment it easily. In contrast, Fenwick tree can be easily coded using array or hash table (unordered_set in STL). Of course, we could also use a regular augmented binary search tree, without the rotation operations to keep the tree balanced, and rely on the fact that if the keys of the input for constructing the binary search tree is random, then the tree height has expectation .