Lazy dynamic programming with arrays in Haskell

In Haskell on November 4, 2011 by Matt Giuca

Dynamic programming is a weirdly named way to speed up (complexity-wise) recursive algorithms. If you are writing a recursive algorithm that ends up calling itself with the same arguments over and over, dynamic programming can solve the problem by “memoizing” the results in a table. Basically, it just means build an array big enough to hold all the possible inputs, then fill in the array as you go. If you ever need the answer to a sub-problem you have already solved, just look it up in the array. The poster-boy for dynamic programming is the knapsack problem: imagine you are a thief who has broken into a jewellery store. You want to steal as much as you can to maximise the total value of the items you’ve stolen, but there is a maximum weight you can carry in your knapsack. Which items do you steal?

The problem with doing this in Haskell is that you can’t modify arrays (well, not easily anyway). Once you create an array, all of its values are fixed. This means you can’t just build an array full of dummy values and go about inserting the actual values later. But, as with many things in Haskell, laziness can save us.

Naïve solution

Let’s start with a simple solution to the unbounded knapsack problem (you can steal as many copies of each item as you like), taken from the Wikipedia article:

-- knapsack [(v0, w0), ..., (vn, wn)] W
knapsack items wmax = m wmax
    m 0 = 0
    m w = maximum $ 0:[vi + m (w - wi) | (vi, wi) <- items, wi <= w]

You call this by passing a list of (v, w) pairs (the value and weight of each item), and the maximum weight the bag can hold, W. You can use the example given in the image at the top of the Wikipedia page:

knapsack [(4, 12), (2, 2), (2, 1), (1, 1), (10, 4)] 15

Which gives the correct answer, 36. It will recursively compute the m function, which tells us the maximum value we can fit into a bag that can carry weight w. However, it is insanely slow because it recomputes sub-problems many times each. To be more formal, it has worst-case time of O(nW), because it is going to test n possibilities to a depth of W. Even this simple example with a maximum weight of 15 takes about four seconds to compute on my laptop. I just bought a new computer yesterday and I’m waiting for it to arrive, which should cut the time down a bit, but in the mean time, we should optimise the algorithm 🙂

Dynamic programming with lazy arrays

It seems like it will be impossible to use dynamic programming in Haskell, because arrays cannot be mutated. But if you think about it, dynamic programming doesn’t actually require that we mutate any values — it merely requires that we have some values that are uncomputed which later become computed. That’s exactly what laziness is all about.

So here’s the plan: we build an array whose indices and elements are the inputs and outputs of the function m. Because of laziness, the array won’t actually contain the outputs of the function — it will just contain thunks which are willing to compute the output if we ask them to. When we recurse, we won’t call a function, we’ll instead pull a value out of the array. If that particular array element has never been requested before, this will cause it to be computed (which might trigger more recursion). If that particular array element has already been found, this will simply return the existing known value. Here we go:

knapsack items wmax = arr ! wmax
    arr = array (0, wmax) [(w, m w) | w <- [0..wmax]]
    m 0 = 0
    m w = maximum $ 0:[vi + arr ! (w - wi) | (vi, wi) <- items, wi <= w]

We first build the array arr as an array of thunks. If any element of the array is requested for the first time, it will call the function m. Then, we try to pull out element wmax, and since it hasn’t been requested, this calls m wmax, which in turn, starts pulling more values out of arr. You can see that this will only ever call m a maximum of wmax times, and furthermore, it won’t necessarily call m for all possible values of w. It now has a worst-case time of O(nW), since each recursion makes n tests, and it will recurse a maximum number of W times. Indeed, it is now very fast in the 15 case. If we change the 15 to, say, 100000, it takes just a few seconds.

Note that all we had to change (besides adding the definition of arr) was replacing all calls to m with a lookup on arr. We should be able to generalise this solution.

A general dynamic programming solution

We would like a solution that allows people to write completely naïve code, such as our first attempt, and automatically turn it into a dynamic programming solution. This is possible in Mercury (see Tabled evaluation) as a special language feature which allows you to say “the runtime should memoize all past outputs of this function”. But we can’t do it without a special language feature (such as in Haskell), because the naïve version recurses on itself. We need to change the recursion slightly to allow our generic tabling code to insert itself in between the recursive calls. So let us go back to our basic solution, but modify m to accept a function argument called “self“, and call that when it recurses. (This is a bit like the Y combinator approach to recursion in Lambda calculus.)

    m _ 0 = 0
    m self w = maximum $ 0:[vi + self (w - wi) | (vi, wi) <- items, wi <= w]

Note that m calls self on the recursive call instead of m. Now it is possible for whatever calls m in the first place to insert some additional code in self. But for now, let’s just make it work as simply as possible (without any tabling). We write a function called table_null:

table_null :: (Ix a) => (a, a) -> ((a -> b) -> a -> b) -> (a -> b)
table_null _ f = let g x = f g x in g

You can pass m to table_null and it will give you back a simple function that computes knapsack values. The first argument to table_null is supposed to represent the domain of the function, but we will ignore it for now. Note that if you pass m as the second argument to table_null, it returns our original m function — when m calls self, table_null‘s g will just call m again. So, this is equivalent to our original, naïve (very slow) knapsack solution:

knapsack items wmax = table_null (0, wmax) m wmax
    m _ 0 = 0
    m self w = maximum $ 0:[vi + self (w - wi) | (vi, wi) <- items, wi <= w]

So what is the point of this? Now we have abstracted away what it means to do recursion. When we use table_null, it just does plain ordinary recursion. But we could make any other table function we like, which can change the way recursion is handled. For example, we could build a table_trace that shows us just how inefficient this approach is by printing out the input to every recursive call:

table_trace :: (Ix a, Show a) => (a, a) -> ((a -> b) -> a -> b) -> (a -> b)
table_trace _ f = let g x = trace (show x) (f g x) in g

If you have knapsack call table_trace, you will see it printing out the values, because table_trace‘s g prints out its input before calling the supplied function.

And now, a general solution to dynamic programming:

table_array :: (Ix a) => (a, a) -> ((a -> b) -> a -> b) -> (a -> b)
  table_array dom f = g
    arr = array dom [(x, f g x) | x <- range dom]
    g x = arr ! x

If you have knapsack call table_array, it will suddenly become efficient. That’s because the g in table_array doesn’t always call f. First, table_array builds an array of thunks as we did above. The thunks inside the array call f the first time they are required. The g function merely picks values out of the array. So anything you use table_array on will automatically memoize the inputs to the function. The only caveat is that you have to know in advance the domain of inputs the function is likely to accept, and it has to be contiguous.


6 Responses to “Lazy dynamic programming with arrays in Haskell”

  1. Nice discussion; nice solution.

    The only limitation is that you need to build the table of thunks before you start. What if the function you want to table takes 5 arguments each covering a ranging of 1000 values? What if it takes an argument with an infinite range (if the table were a list, you could handle it, but I’ve never heard of an infinite array in Haskell)? How could you use this to compute fib?

    • So if you had some really complex or wide-ranging inputs, and you were planning to sparsely look up a relatively small percentage of the ones you need, then clearly you want more of a hash or tree like structure than an array.

      I thought about adding a suggestion like that to the end, but then I realised that it wouldn’t work with this approach — precisely as you say, you need to build all the thunks before you begin. I suppose all of the dynamic programming that uses an array does need to be relatively dense because it does make the problem linear in the size of the entire array, in both space and time (e.g., knapsack, while it won’t necessarily query all values, it will query every multiple of the lowest common multiple of the weights, so if you take wiki’s advice and divide them all by the LCM first, then it will completely cover the array).

      Couldn’t you compute fib with this? You can’t make an infinite array, but if you wanted to compute fib(N), you could make an array of size N. Interestingly, you could make an infinite list of thunks, but then it would become quadratic as the lookup time would be linear.

      • Yes, you could implement fib this way by initialising the array for each call. But then it wouldn’t remember results between non-recursive calls to fib. The ideal data structure for this would be a lazy infinite array that you allocate by specifying only the initial value for each array element as a function of its index. When you access any element for the first time, it materialises the value and stores it, so thereafter when you access that element, it only requires a lookup.

  2. Yeah, that’s the same problem for knapsack, above. But there isn’t such a thing as an “infinite array”, though I suppose if your array data type was more like a vector which expanded itself as necessary, that would work. That’s an interesting idea for a data structure (which would need to be implemented primitively) — an expandable potentially infinite array with lazy elements. It can’t be written to, only read. But it allocates finite memory based on the first element read, and then whenever you read it, if you read further than ever before, it reallocates and expands the array.

    Anyway, point is that using arrays in C gives the same limitation as the Haskell approach (and the same space and time complexity). So it isn’t fully lazy as you might expect in Haskell, but at least it’s not worse.

  3. […] came across a cool blog post by Matt Giuca where he suggests that we need to think of things which have been computed and those which […]

  4. actually, dynamic programming is a little bit more than just memoizing all the function computations. It generally also involves analyzing the *structure* of the recursion in order to find a minimal set of values which need to be reused. This information is then used to formulate the computation as an accumulation in the reversed data dependency direction from the recursive formulation. This avoids the problems pointed out by the commenters concerning over-initialization (which will kill your complexity guarantees, just as surely as the repeated computations in the recursive formulation).

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: