akrish

I’m continuing to learn Haskell in my free time and as I like to learn by doing, I’ve been implementing common algorithms in the language. Last week I tackled depth first search and found it to be surprisingly challenging and consequently really educational. Today, I decided to try my hand at dynamic programming. I chose dynamic programming because these algorithms hinge on memoization, and I wanted to figure out how to do that in a purely functional way. Also, I thought that I could do some performance analysis, as implementing recursive, un-memoized algorithms to solve these problems seemed to be quite easy.

The first algorithm I tried to implement was the standard edit distance algorithm. It’s a classic DP algorithm and although no longer used much in practice, it is the basis for several modern techniques used in sequence alignment, so I figured it would be an interesting exercise. I spent an hour or so thinking about how to approach this and even started programming a solution, but I quickly got bogged down in maintaining a matrix of subproblem solutions. I found implementing all of these matrix operations to be quite mundane and decided to start with something a bit less tedious.

So I chose another classic DP problem, Longest Increasing Subsequence. Briefly, the problem formulation is: Given a list of (lets say) positive, unique integers, find the length of the longest subsequence such that the subsequence is strictly increasing. For example, given [1, 4, 2, 3, 9], the longest increasing subsequence would be [1, 2, 3, 9]. Read the wikipedia article for more details. I chose this problem because you don’t have to maintain a matrix of subproblem solutions, but rather a list, and the problem is quite a bit simpler, allowing me to focus on the algorithmic elements.

So, I started out by implementing a naive un-memoized algorithm, then went on to implement the dynamic programming algorithm that’s based on the one on the wikipedia page, but in Haskell it doesn’t really look the same at all. I implemented both so that I could do run some performance metrics which I’ll discuss about after presenting my code. First some auxiliary stuff:

[code lang="lisp"]
-- return the item in the list indexed by i --
listGet :: [Int] -> Int -> Int
listGet lst i =
if i == 0
then (head lst)
else
listGet (tail lst) (i-1)

-- return the index of item i in the list --
getIndex :: [Int] -> Int -> Int
getIndex lst i =
if (head lst) == i
then 0
else
1 + (getIndex (tail lst) i)

-- given a list of integers, return the largest value in the list --
listMax :: [Int] -> Int
listMax [] = 0
listMax lst =
if (length lst)==1
then (head lst)
else max (head lst) (listMax (tail lst))
[/code]

I admit I should have used arrays for this but I kind of started using lists before having to implement these and was too lazy to change it. I would change it now but I’m worried about breaking things right before posting and since I’m probably (read definitely) never going to look at this again it really doesn’t matter to me.

Anyway, so we have three functions that are all pretty straightforward. One essentially treats a list as an array and gets the i-th element from the array. Another does the reverse, it gets the index of an item in the array. Finally the third computes the max over a list.

My first algorithm is a naive approach that is just recursive with no memoization:

[code lang="lisp"]
-- compute the longest increasing subsequence starting from each index in the list --
longestIncreasingSubsequence :: [Int] -> Int -> [Int]
longestIncreasingSubsequence [] _ = []
longestIncreasingSubsequence lst i =
if (length lst) == i+1
then [1]
else
[(1+(listMax
(map (\x -> head (longestIncreasingSubsequence lst (getIndex lst x)))
(filter (\x -> (listGet lst i) < x)
(filter (\x -> (getIndex lst x) > i) lst)))))]
++ (longestIncreasingSubsequence lst (i+1))

-- compute the longest increasing subsequence of the list --
lis2 :: [Int] -> Int
lis2 lst =
listMax (longestIncreasingSubsequence lst 0)
[/code]

Ok so longestIncreasingSubsequence is a bit confusing and the comment is kind of inaccurate. If you call it with i == 0, then it’ll return a list where each j-th element is the value of the longest increasing subsequence that starts from and contains the jth item of lst. It basically computes the solution to all of the subproblems, but it doesn’t compute them by memoization. It’s similar if you start with i != 0 but the list only has a suffix of the subproblems. It computes these results by following the standard recursive formula: look at each j>i such that l[i]

This was a good exercise but not really what I wanted to do. So I tested that a bit and made sure it was correct, and proceeded to think about the dynamic programming algorithm. I found it a bit difficult to actually memoize the solutions, but here's what I came up with:

[code lang="lisp"]
-- given an array of subproblem, find the largest increasing subsequence starting at i --
dplis :: [Int] -> Int -> Array Int Int -> Int
dplis lst int arr =
listMax
((map (\x -> 1 + arr ! (getIndex lst x))
(filter (\x -> (listGet lst int) < x)
(filter (\x -> (getIndex lst x) > int) lst))) ++ [1])

– compute largest increasing subsequence starting at int –
getDPScores :: [Int] -> Int -> Array Int Int
getDPScores lst int =
if (length lst) == int+1
then array (0, ((length lst)-1)) ([(i, 0) | i <- [0..((length lst)-2)]]++[(((length lst)-1), 1)])
else
let arr = getDPScores lst (int+1) in
arr // [(int, dplis lst int arr)]

– compute largest increasing subsequence of list –
lis :: [Int] -> Int
lis lst =
listMax (elems (getDPScores lst 0))
[/code]

So I hope the functions make sense. dplis just does lookups in the array parameter to compute the longest increasing subsequence starting at and including i. getDPScores builds the array of subproblem solutions, starting from the smallest subproblems and relying on dplis. Finally lis is just a wrapper for calling getDPScores with the right parameters.

One noteworthy thing is the flow of the program. It's tail recursive. So the first thing we do is get to the bottom of the recursion of getDPScores, then we really do our computation as we work our way out of that recursion stack. Each call to getDPScores after returning from recursion still has to call dplis. Of course there are a lot of tail recursive algorithm, but I didn't really expect dynamic programming to naturally flow like this. Typically dynamic programming is done iteratively, but I suppose we can't do that here.

So those are my two implementation, now on to performance stuff. I decided to experiment with this just to see how amazing dynamic programming really is. If you notice the first algorithm I implemented is definitely exponential (in the worst case, imagine a ascending list of length n, then each call has to spawn O(n) smaller calls). I think it's O(n^n) but I haven't thought too much about it so please correct me if I'm wrong. Anyway, the second algorithm is O(n^2), you just recurse straight down and pop all the way back up. At each recursive call you may have to look at O(n) neighbors and there are only n recursive calls.

Given the asymptotic running times, I expect the second algorithm to be far superior to the first in actual running time. Just for fun I decided to run this experiment and plot my results. I generated random lists of sizes ranging from 5 to 70 in increments of 5 and ran each algorithm on the same lists. Here are the results:

They are pretty much what I expected. Notice that the y axis is log-scaled, so that for input length 70, the exponential algorithm took around 10 minutes (600 seconds) where the dynamic programming one takes well under a second. Really amazing how dynamic programming helps that much. I had originally intended to run on larger input sizes but I simply couldn't as the exponential algorithm wouldn't terminate.

As I did in my last Haskell post, I did a quick search for other people who had looked at dynamic programming in Haskell. I found this article which walks through an implementation of the knapsack algorithm and it appears to be much more elegant than mine. Of course you should remember that this I'm still a Haskell newbie so I don't have mastery over most of the features and I really don't expect myself to be able to implement elegant solutions. However, he does mention that dynamic programming is quite natural in functional languages and that you can take advantage of lazy evaluation, which I can understand. Also the knapsack algorithm uses a matrix of subproblems, so it's slightly different and more interesting.

I also found this article that implements the Needleman/Wunsch algorithm for global sequence alignment, another classic dynamic programming algorithm. This was also an interesting read, and he used some language features that I haven't used and barely knew existed so it was definitely educational. Also, this algorithm also uses a matrix of subproblem solutions so I'm thinking that my next experiment will be with more one of these dynamic programming algorithms.

So, In conclusion, I'm feeling more confident in my functional thinking and my Haskell programming. In my previous post, Jedai commented that I shouldn't consider arrays as cheating, and so I have adopted them into my repertoire here. Further, as proof of the importance of algorithm design and asymptotic efficiency, I generated some performance numbers to show how essential it is to have efficient algorithms.

Anyway... there's a leak in my ceiling and water is dripping all over the place so, until next time!

Tags: haskell cs programming dp

This entry was posted on Sunday, June 28th, 2009 at 2:01 pm and is filed under cs. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.