Dynamic Programming in Haskell

Submitted by mrd on Thu, 02/22/2007 - 8:57pm.

Dynamic Programming is an algorithm design strategy which can be essentially described as breaking a problem down into subproblems and re-using previously computed results. Sometimes, a problem which appears to take exponential time to solve can be reformulated using a Dynamic Programming approach in which it becomes tractable. The benefit is especially clear when the subproblem solutions overlap considerably. The technique of memoization is a major time-saver in those cases.

A common Haskell newbie question is: how do I memoize? At first, it appears to be a very difficult problem, because access to mutable arrays and hashtables is restricted. It is important to realize that lazy evaluation is actually memoization itself and can be leveraged in that way for the purposes of Dynamic Programming. In fact, as a result, the expression of these algorithms can be more natural and lucid in Haskell than in a strict language.

Here, I am going to examine the classic ``knapsack problem.'' Given a number of items, their values, and their sizes -- what is the best combination of items that you can fit in a limited-size knapsack?

> module Knapsack where > import Control.Monad > import Data.Array > import Data.List > import Test.QuickCheck
I am going to represent items with this data type. Essentially, it is just a tuple of the item itself, its value, and its size. > data Item a = Item { item :: a, > itemValue :: Int, > itemSize :: Int } > deriving (Eq, Show, Ord)
Cells will be used both to represent the solution to the knapsack problem, and as individual cells in the matrix for the Dynamic Programming algorithm. It is a pair consisting of: summed values, and the items in the sack. > data Cell a = Cell (Int, [Item a]) > deriving (Eq, Show, Ord)
This powerset definition is a very neat use of the List monad (which I pulled from the Haskell wiki). You can think of it as saying: for each element in the list, half the possible subsets will include it, and half will not. > powerset :: [a] -> [[a]] > powerset = filterM (const [True, False])
brutepack considers the powerset of the items, cutting out those subsets which are too large in size, and picking the most valuable subset left. As you might figure, this is going to run in O(2^n) thanks to the use of powerset. The definition should be simple to understand and will provide a sound basis for testing the Dynamic Programming alternative. > brutepack :: (Ord a) => Int -> [Item a] -> Cell a > brutepack size items = maximum [ > cellOf subset | > subset <- powerset items, itemsFit subset > ] > where > itemsFit items = sum (map itemSize items) <= size > cellOf items = Cell (sum (map itemValue items), items)

The Dynamic Programming algorithm is as follows:

Consider a matrix where the rows are indexed by size and the columns by items. The rows range from 1 to the size of the knapsack, and the columns are one-to-one with the items in the list.

value = $30     $20    $40
 size =  4       3      5
     4|        v(2,4)<--.
      |                 |
     5|                 |
      |                 |
     6|                 |
      |                 |
     7|                 |
      |              <--|    
     8|        v(2,8) v(3,8)

This is a diagram where the knapsack has a maximum size allowance of 8, and we want to stuff some animals in it. Each element in the matrix is going to tell us the best value of items by size. That means the answer to the whole problem is going to be found in v(3,8) which is the bottom-rightmost corner.

The value of any one cell in the matrix will be decided by whether it is worthwhile to add that item to the sack or not. In the v(3,8) cell it compares the v(2,8) cell to the left to the v(2,4) cell up above. The v(2,8) cell has no room for the bear, and the v(2,4) cell represents the situation where the bear will fit. So the question, after determining if the bear will fit in the bag at all, is: is value of bear + v(2,4) better than v(2,8)?

This definition of v lends itself to a very nice recursive formulation.

              / v(m-1,n)                    if (s_n) > n
    v(m,n) = -|      / v(m-1,n)
              /      \
              \ max -|                      otherwise
                     \ v(m-1,n-(s_n)) + v_n
where s_n is the size of item n and v_n is its value.

A typical implementation of this algorithm might initialize a 2D array to some value representing ``undefined.'' In Haskell, we can initialize the entire array to the correct value directly and recursively because it will not be computed until needed. All that is necessary is to express the data-dependencies, and the order of evaluation will take care of itself.

This code takes the algorithm a little further by tracking a field in each cell that contains a list of the items in the sack at that point.

> dynapack :: Int -> [Item a] -> Cell a > dynapack size items = valOf noItems size > where > noItems = length items > itemsArr = listArray (1,noItems) items > itemNo n = itemsArr ! n > > table = array ((1,1),(noItems,size)) $ > [ > ((m, n), cell m n) | > m <- [1..noItems], > n <- [1..size] > ] > > valOf m n > | m < 1 || n < 1 = Cell (0, []) > | otherwise = table ! (m,n) > > cell m n = > case itemNo m of > i@(Item _ v s) > | s > n || > vL >= vU + v -> Cell (vL , isL) > | otherwise -> Cell (vU + v, i:isU) > where > Cell (vL, isL) = valOf (m - 1) n > Cell (vU, isU) = valOf (m - 1) (n - s)

The matrix is defined in the variable table and valOf is our function v here. This definition very naturally follows from the algorithmic description because there is no problem with self-reference when defining cells in the array. In a strict language, the programmer would have to manually check for the presence of values and fill in the table.


It's important to be confident that the algorithm is coded correctly. Let's see if brutepack and dynapack agree on test inputs. I will define my own simple data type to customize the randomly generated test data.

> newtype TestItems = TestItems [(Int, Int, Int)] > deriving (Eq, Show, Ord)
> nubWithKey k = nubBy (\a b -> k a == k b) > fst3 (a,b,c) = a > tripleToItem (i,v,s) = Item i v s

The Arbitrary class expects you to define your customly generated test data in the Gen monad. QuickCheck provides a number of handy combinators, and of course you can use normal monadic functions.

sized is a QuickCheck combinator which binds the generator's notion of "size" to a parameter of the supplied function. This notion of "size" is a hint to test-data creators that QuickCheck wants data on the "order of" that size. Of course, what "size" means can be freely interpreted by the author of the function, in this case I am using it for a couple purposes.

The basic idea is simply: create a list of randomly generated tuples of length "size", and choose values and item-sizes randomly from (1, "size"). Notice how the randomly generated tuple is replicated with the monadic combinator replicateM. Then, before returning, just make sure that there are no "repeated items" by running nubWithKey fst3 over the generated list. That will cut out any items with the same name as previous items.

> instance Arbitrary TestItems where > arbitrary = sized $ \n -> do > items <- replicateM n > $ do > i <- arbitrary > v <- choose (1, n) > s <- choose (1, n) > return $ (i, v, s) > return . TestItems . nubWithKey fst3 $ items

With an Arbitrary instance, we can now define a property: it extracts the tuples and creates Items out of them, then tries out both algorithms for equivalence. Note that I am only checking the final values, not the actual items, because there may be more than one solution of the same value.

> prop_effectivePacking (TestItems items) = v1 == v2 > where items' = map tripleToItem items > Cell (v1,_) = brutepack 16 items' > Cell (v2,_) = dynapack 16 items'
Knapsack> verboseCheck prop_effectivePacking 0: TestItems [(1,2,3),(3,2,2)] 1: TestItems [] 2: TestItems [(1,1,1)] 3: TestItems [(2,1,2),(-1,2,3),(0,2,3)] 4: TestItems [(-2,2,2),(-1,1,1),(3,2,3),(4,2,2)] ...

It will progressively check larger "size" samples and you will notice that the brute force algorithm is going to start dragging down performance mightily. On my computer, in the ghci interpreter, brutepack on even just 20 items takes 10 seconds; while dynapack takes almost no time at all.


Algorithms making use of Dynamic Programming techniques are often expressed in the literature in an imperative style. I have demonstrated an example of one such algorithm in a functional language without resorting to any imperative features. The result is a natural recursive expression of the algorithm, that has its advantage from the use of lazy evaluation.

Submitted by Cale Gibbard on Mon, 02/26/2007 - 12:06pm.

Something seems to have gone wrong with the formatting here.

Submitted by mrd on Mon, 02/26/2007 - 3:09pm.
How so?
Submitted by mrd on Mon, 02/26/2007 - 6:53pm.

Note: bhb on reddit pointed out a mistake in the recursive formula diagram, now corrected.

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.