Week 9: Notes

2-dimensional dynamic programming

The dynamic programming problems we saw last week had a 1-dimensional structure. That means that each of those problems could be solved by writing a recursive function with 1 argument. The naive recursive solution was exponentially inefficient, since it solved the same subproblems over and over again. We found that we could solve these problems much more efficiently by filling in a 1-dimensional array in a certain order. This array recorded the solution to each subproblem, so that it was immediately available to use in solving other, larger, subproblems.

For example, the naive recursive solution to the rod-cutting problem looked like this:

  // Return the best price for cutting a rod of length n, given a table
  // with the prices for pieces from lengths 1 .. n.
  int profit(int n, int[] prices) {
      int best = 0;
      for (int i = 1 ; i <= n ; ++i)
          best = Max(best, prices[i] + profit(n - i, prices));
      return best;
  }

Notice that profit() is a function of a single argument "int n" (ignoring the constant array prices[], which could just be stored as a field outside the function). To compute the result more efficiently, we used a bottom-up approach that filled in a 1-dimensional array called best[], which held the best possible profit for a rod of size i, i.e. the value that would be returned by profit(i, prices).

We will now study slightly more challenging dynamic programming problems that have a 2-dimensional structure. For these problems, a naive recursive solution will be a function with 2 arguments. Once again, this naive solution will run in exponential time. We can solve these problems more efficiently using bottom-up dynamic programming in which we fill in a 2-dimensional array or array-like data structure. We will need to be careful about the order in which we fill in the array elements. The solution to any subproblem instance will depend on the solutions to smaller instances, so we need to fill in the array in some order that guarantees that the dependent solutions will already be available when we compute any subproblem solution.

number of ways to form a sum

As a first example, recall that last week we saw how to write a recursive function that can print out all ways to form a sum from Czech coins. We wrote two versions of the function. In the first version, order matters, so e.g. 1 + 2 and 2 + 1 are considered to be distinct sums. In the second version, order does not matter, so e.g. 1 + 2 and 2 + 1 are the same sum, and we only want to print one of these. We accomplished that by forcing each sum to contain a non-increasing series of values. See last week's lecture notes for the programs we wrote.

Now consider the related problem of computing only the number of possible ways to form a sum. Let's first assume that order matters, so we want to count all permutations of each sum (e.g. 1 + 2 and 2 + 1) separately. Here's a recursive function to perform this computation:

int[] coins = { 1, 2, 5, 10, 20, 50 };

long ways(int n) {
    if (n < 0)
        return 0;
    else if (n == 0)
        return 1;   // 1 way to sum to 0: the empty set
    else {
        long w = 0;
        foreach (int c in coins)
            w += ways(n - c);
        return w;
    }
}

Let's run this function for an increasing sequence of values of n:

n = 0: ways = 1
n = 1: ways = 1
n = 2: ways = 2
n = 3: ways = 3
n = 4: ways = 5
n = 5: ways = 9
...

n = 30: ways = 5,878,212
n = 31: ways = 10,050,981
n = 32: ways = 17,185,883
n = 33: ways = 29,385,638
n = 34: ways = 50,245,647
...

We see that the function's value increases exponentially. Also, it becomes very slow for values of n that are more than about 30.

The function is slow because it's solving the same subproblems over and over again. We can make it much more efficient using dynamic programming. This is a 1-dimensional dynamic programming problem, since ways() takes a single argument. Here is a dynamic programming solution:

int[] coins = { 1, 2, 5, 10, 20, 50 };

long ways_dp(int n) {
    long[] ways = new long[n + 1];

    ways[0] = 1;
    for (int k = 1 ; k <= n ; ++k) {
        long w = 0;
        foreach (int c in coins)
            if (k - c >= 0)
                w += ways[k - c];
        ways[k] = w;
    }

    return ways[n];
}

This function will be far more efficient than the previous version. In fact ways_dp(N) will run in O(N), since the outer for loop iterates N times and the inner foreach loop iterates only a constant number of times.

Now consider computing the number of ways to form a sum when order does not matter: 1 + 2 and 2 + 1 are the same sum, so we only want to count one of these. We can compute this number recursively. As with the function we wrote last week, we'll add an extra argument 'max' that allows us to limit the value of any coin we can choose from that point onward. Here is the recursive function:

int[] coins = { 1, 2, 5, 10, 20, 50 };

long _ways2(int n, int max) {
    if (n < 0)
        return 0;
    else if (n == 0)
        return 1;
    else {
        long w = 0;
        foreach (int c in coins)
            if (c <= max)
                w += _ways2(n - c, c);
        return w;
    }
}

long ways2(int n) => _ways2(n, 50);

Let's run it for an increasing sequence of values of n:

n = 0: ways2 = 1
n = 1: ways2 = 1
n = 2: ways2 = 2
n = 3: ways2 = 2
n = 4: ways2 = 3
n = 5: ways2 = 4
n = 6: ways2 = 5
n = 7: ways2 = 6
...

n = 300: ways2 = 393,119
n = 301: ways2 = 398,610
n = 302: ways2 = 404,781
n = 303: ways2 = 410,272
n = 304: ways2 = 416,443

...

This function increases much more slowly than the previous function, since there are many fewer ways to form a given sum when all permutations of a sum are considered to be the same. However this function will also become very slow for large values of n, so we'd like to make it more efficient using dynamic programming.

The recursive function _ways2() takes two arguments, indicating that this will be a 2-dimensional dynamic programming problem. We want to build a 2-dimensional array that holds the value of _ways(n, max) for any possible values of n and max.

The largest possible value of max is 50, the highest coin value. However it would be wasteful to have an array of size N x 50, since max will always be one of the 6 possible coin values. So let's first rewrite _ways2() so that its second argument max_index is the index of one of the values in the coins[] array. Then max_index will be an integer in the range 0 .. 5, since there are 6 possible coin values. Here is the rewritten recursive function:

int[] coins = { 1, 2, 5, 10, 20, 50 };

long _ways2(int n, int max_index) {
    if (n < 0)
        return 0;
    else if (n == 0)
        return 1;
    else {
        long w = 0;
        for (int i = 0 ; i <= max_index ; ++i)
            w += _ways2(n - coins[i], i);
        return w;
    }
}

long ways2(int n) => _ways2(n, coins.Length - 1);

Now we can transform this function into an efficient solution using dynamic programming. Here is the code:

long ways2_dp(int n) {
    // ways[k, i] is the number of ways to form the sum k using coins
    // whose index is at most i.
    long[,] ways = new long[n + 1, coins.Length];

    for (int i = 0 ; i < coins.Length ; ++i)
        ways[0, i] = 1;    // 1 way to make sum of 0

    for (int k = 1 ; k <= n ; ++k)
        for (int max_index = 0 ; max_index <= coins.Length - 1 ; ++max_index) {
            long w = 0;
            for (int i = 0 ; i <= max_index ; ++i)
                if (k - coins[i] >= 0)
                    w += ways[k - coins[i], i];
            ways[k, max_index] = w;
        }

    return ways[n, coins.Length - 1];
}

Notice the triple nested for loop above. The outer loops 'for (int k = 1 ...' and 'for (int max = 0 ...' iterate over all cells of our 2-dimensional table. The inner loop 'for (int i = 0 ; i <= max_index ...' is performing the same calculation that we made in the recursive function _ways2(). This is a 2-dimensional solution but will still run in O(N), since the table width is constant and since only one of the nested for loops iterates N times.

In any 2-dimensional dynamic programming problem we need to be careful about the order in which we fill in the table cells, to be sure that each cell's value will be computed before we first need it. In this problem, the calculation of ways[k, max_index] uses the values of table cells ways[k - coins[i], i] for various values of i. Notice that k - coins[i] < k, and i <= max_index. So we can fill in the array one row at a time, using a double loop in which k and max_index both increase, and we'll have each cell value by the time we need it. As we'll see in our next dynamic programming problem below, sometimes we'll need to use a different order.

longest palindromic subsequence

Consider the problem of finding the longest palindromic subsequence (LPS) of a given string. For example, consider the string "WATERMELON STEW". The string "WTERETW" is a longest palindromic subsequence:

W A T E R M E L O N   S T E W
W   T E R   E           T   W

That's because this subsequence is a palindrome (i.e. it reads the same forwards and backwards), and there is no longer palindromic subsequence in the string. Note that the longest palindromic subsequence is not necessarily unique: "WTEMETW" is another possibility.

We'll first write a recursive solution to this problem, then transform it to a more efficient solution using 2-dimensional dynamic programming.

In our discussion it will be convenient to use the notation s[i ... j] to denote the substring of s consisting of all characters of s from s[i] up to and including s[j]. Note the triple dots in this notation. This is distinct from the C# expression s[i .. j] (written with only two dots), which is a substring that does not include s[j].

Now suppose that the input string s has length N, with character indices from 0 to (N - 1). We want to construct a recursive function P(i, j) that for any i and j computes some longest palindromic subsequence of s[i ... j]. Once we have P, we can simply call P(0, s.Length - 1) to find a longest palindromic subsequence of the entire string s.

Our base case is straightforward: if i = j, then we can choose P(i, j) = s[i], since any 1-character string is a palindrome. In fact another base case is even smaller: if j < i, then s[i ... j] is empty, so P(i, j) = "", i.e. the empty string. In either of these cases P(i, j) = s[i ... j], which is the same as s[i .. (j + 1)] in C#.

Now suppose that i < j, which will be the recursive case. We may consider two possibilities:

Since we have found a recursive pattern, we can now implement the function P, which we will now call _lps:

string longer(string s, string t) =>
    s.Length > t.Length ? s : t;

// compute lps of s[i ... j ]
string _lps(string s, int i, int j) {
    if (j <= i) // s[i ... j] has length 0 or 1
        return s[i .. (j + 1)];
    
    if (s[i] == s[j])
        return s[i] + _lps(s, i + 1, j - 1) + s[j];
    else
        return longer(_lps(s, i, j - 1), _lps(s, i + 1, j));
}

string lps(string s) => _lps(s, 0, s.Length - 1);

It works:

 WriteLine(lps("WATERMELONSTEW"));  // writes WTEMETW

However this function will run in exponential time.

To solve this problem more efficiently, we will use bottom-up dynamic programming to fill in a 2-dimensional array lps[], where lps[i, j] = _lps(i, j).

Let's think a bit about the structure of the table we are building. Assume as usual that i is a row number and j is a column number. If i = j, then lps[i, j] is on the diagonal, and will contain a single-character string holding just the letter s[i]. That corresponds to the first base case we described above. If i > j, then lps[i, j] is below the diagonal and will be the empty string, since s[i ... j] is empty. We don't even need to fill in those table cells. The interesting cells are above the diagonal, where i < j and s[i ... j] has at least two characters.

We need to be careful about the order in which we fill in the array. Specifically, when we compute lps[i, j] we must already know lps[i + 1, j] and lps[i, j – 1] and lps[i + 1, j – 1]. If we fill in rows from top to bottom, we will have a problem, since at the time we compute lps[i, j] we will not yet know lps[i + 1, j] or lps[i + 1, j - 1]. So we will fill in rows from bottom to top. We can fill in each row from left to right, since when we compute lps[i, j] we don't need the values of any table cells with higher values of j.

As we fill in each row i, we'll begin with the diagonal element lps[i, i], and proceed rightward from there. Here is the code:

string longer(string s, string t) =>
    s.Length > t.Length ? s : t;

string lps_dp(string s) {
    string[,] lps = new string[s.Length, s.Length];

    for (int i = s.Length - 1 ; i >= 0 ; --i) {
        lps[i, i] = "" + s[i];
        for (int j = i + 1 ; j < s.Length ; ++j) {
            if (s[i] == s[j])
                lps[i, j] = s[i] + lps[i + 1, j - 1] + s[j];
            else
                lps[i, j] = longer(lps[i, j - 1], lps[i + 1, j]);
        }
    }

    return lps[0, s.Length - 1];
}

How efficient is this function? Consider running lps_dp(S), where S is a string of length N. In the worst case, S itself is a palindrome, and then the longest palindromic subsequence of S is S itself. Then we're filling in the upper half of a table of N x N elements, and each table element will be a string whose length may have as many as N characters. The worst-case time will be O(N³).

We could make this function more efficient. As one possibility, instead of storing a string in each table cell lps[i, j], we could store an integer which is the length of the longest palindromic subsequence of s[i ... j]. Then, after we've built the table we could have a loop that iterates over the table and reconstructs a longest palindromic subsequence. For details, see last year's lecture notes, where I included code for this sort of solution.

As another possibility, we could make the function more efficient by storing strings in a way that allows us to prepend or append characters in O(1). This is possible using linked lists of characters, and can also lead to a solution that runs in O(N²). (This sort of solution will be easier in a functional language such as Haskell, which you'll learn if you take Non-Procedural Programming next year.)