Lambda the Ultimate

I think that this short paper by Peter Norvig is a very nice demonstration of dynamic programming via memoization in Lisp: http://acl.ldc.upenn.edu/J/J91/J91-1004.pdf

Here's a really basic subset I just hacked up in (not-quite-portable) Common Lisp:

(defmacro define-recurrence (name &body equations)
  (let ((LHSs (mapcar #'first equations))
	(RHSs (mapcar #'third equations)))
    (let ((n (car (last LHSs))))
      (expand-recurrence name n LHSs RHSs))))


(defun expand-recurrence (name n LHSs RHSs)
  (labels ((expand-equation (lhs rhs)
	     `(,(expand-lhs lhs) ,(expand-rhs rhs)))
	   (expand-lhs (lhs)
	     `(= ,n ,lhs))
	   (expand-rhs (rhs)
	     (cond ((atom rhs) rhs)
		   ((eq (first rhs) name)
		    `(aref ,@(mapcar #'expand-rhs rhs)))
		   (t (mapcar #'expand-rhs rhs)))))
    `(defun ,name (,n)
      (let ((,name (make-array (list (+ ,n 1)))))
	(dotimes (,n (+ ,n 1))	;FIXME: may assign instead of bind
	  (setf (aref ,name ,n)
		(cond ,@(mapcar #'expand-equation LHSs RHSs))))
	(aref ,name ,n)))))


(define-recurrence fib
  (0 -> 0)
  (1 -> 1)
  (n -> (+ (fib (- n 1))
	   (fib (- n 2)))))

(defmacro defun-recurrence (name (n) &body body)
  `(defun ,name (,n)
     (let ((,name (make-array (+ ,n 1))))
       (loop for ,n from 0 to ,n
             do (setf (aref ,name ,n)
                      (flet ((,name (,n) (aref ,name ,n)))
                        ,@body)))
       (aref ,name ,n))))


(defun-recurrence fib (n)
  (if (< n 2)
      n
      (+ (fib (- n 1))
(fib (- n 2)))))

So if we generalized that to multidimensional arrays, would DPROG be doing anything more? I'm not sure.

At http://www.cs.sunysb.edu/~liu/ there's a paper "Dynamic programming via static incrementalization" about analyzing/transforming code in ordinary languages in an intelligent way -- I only skimmed the first few pages but it looks like it's closer to Ehud's idea of inferring how the DP algorithm should work...

And at http://accesscom.com/~darius/hacks/screenfuls/screen3.html is a fancier version of the Lisp above. It can handle all the examples from the DPROG tutorial -- e.g.

  M[n,k] = select{A[0]                            when n=1,
                  sum{ A[i] where 0 <= i < n }    when k=1,
                  min{max{M[i,k-1], sum{ A[j] where i <= j < n } }
                      where 1 <= i <= n } }
where 1 <= n <= N and 1 <= k <= K

becomes

(defun-recurrence M (n k) (1 1)
  (cond ((= n 1) (A 0))
        ((= k 1) (loop for i from 0 below n 
                       sum (A i)))
        (t (loop for i from 1 to n
                 minimize (max (M i (- k 1))
                               (loop for j from i below n
sum (A j)))))))

It's so much easier to extend the right language than create a new one! Of course, if you start getting into things like fancy analyses to detect which subset of the values you'll need, it may be just as much work in the end, but at least you got something usable quickly.

If I wanted to use this for an algorithm like FIB that only needs the last few values, not a table of all values, I'd modify the TABULATE macro in the code linked above. (Basically turn (AREF a i) into (AREF a (MOD i 2)).) You could make an open implementation with hooks for different changes, like hand-rolled AOP -- I won't bother because then the code would no longer fit on one screen. :-)

Finally, it's possible to do this with a C++ library after all -- code using it would look like

inline int compute_fib(Recurrer<int, int>& fib, int n) {
  return n < 2 ? n : fib(n-1) + fib(n-2);
}


int fib(int n) {
  return recur(compute_fib, n);
}

I see that you are discussing DPROG and I would like to respond to a few comments on it. First of all, the memorization versus dynamic programming choice: The main reason for using dynamic programming and not memorization for my purposes is that I need to extend the dynamic programming algorithms generated with garbage collcection, so to speak, of the tables filled in during the algorithm. The data sets I work with are too large to allocate N*N or N^3 matrices, so I only want to allocate and use the columns and rows that I strictly need for the algorithm (and then use Hirschbergs algorithm to re-generate paths through the matrices when I need them, or just the the final value when I only need that). With dynamic programming this is simple to achieve, with memorization I wouldn't know how to approach it.

This leads me to the second comment: Why write a new compiler instead of using lisp macros or C++ templates? The main reason for writing DPROG is to get the compiler to automatically infer the number of columns and/or rows needed at any given time to complete the algorithm, and provide hooks to Hirschberg's algorithm to simplify using it. Calculating dependencies between matrices and matrix-cells is what I am currently working on. At least, it is what I was working on at the start of this week and is what I will be working on when I return from the conference I am currently attending. This dependency calculation is an important part of the compiler, and will be added before moving from alpha to beta. Using a library in another programming language can be just as easy as expressing the recurrence in DPROG, at least for simpler recurrences, but using the DPROG compiler you can (hopefully) check that the recurrence is well-founded and infer which elements in the algorithm you need to remember and which you can delete as the algorithm progresses.

Yours, /mailund

Sparse arrays are commonly used for this purpose.

With dynamic programming this is simple to achieve, with memoization I wouldn't know how to approach it.

You could use sparse arrays with a solution like the one in Darius Bacon's second post in this thread.

BTW, I'm not arguing against the use of a DSL. However, it seems to me that an embedded DSL could be more useful than a standalone one in a case like this. If this were embedded in Scheme or Lisp, you could easily analyze recurrence definitions. And it would have saved you a lot of work! ;)

P.S. I notice the Norvig paper referenced above uses hash tables, which address the sparseness issue.

Sparse arrays are commonly used for this purpose.

I haven't had time to read the Norvig paper, so what I am about to say might be way off...

As I understand the suggestion of using sparce arrays, you are saying that I can reduce memory usage by storing only the entries in the arrays that I need. That is true, but what if I really need all entries, just not at the same time in the algorithm? Typically, I need to fill in a matrix one row or one column at a time, and when filling a row/column I only refer to the current row/column or the previous row/column. All rows/columns further back are no longer needed, so I want to re-use the memory in those rows/columns.

If, for instance, I fill in the matrix one row at a time, only refering back a single row, I can adress into the matrix as M[i,j] == M + N*(i mod 2) + j (where N is the rowlength).

The problem with memorization (as I understand the sparse arrays to be) is that it saves me from using memory for cells that will never be used, but it does not allow me to recognize when a cell is no longer of use and to re-claim the memory for another cell.

In the applications I am interested in, I typically need all cells but not at the same time, so memorization is not the right choice.

BTW, I'm not arguing against the use of a DSL. However, it seems to me that an embedded DSL could be more useful than a standalone one in a case like this. If this were embedded in Scheme or Lisp, you could easily analyze recurrence definitions. And it would have saved you a lot of work! ;)

I absolutely agree that embedding the DSL in an existing language is better than what I do with DPROG, unfortunatly I need to embed it in C++ and I couldn't think of a way to do this easily--especially if I wanted the analysis of the program to automatically determine the order to fill in cells and to reduce memory usage by reclaiming memory for storing cells.

You are probably right about it being easier in Lisp. My lisp experience is limited to hacking a bit of elisp and playing around with guile a few years ago, so I wouldn't know how to approach it, though.

But as I see it, the only thing I would really save in effort by programming this in lisp is the parsing and code generation (since the input program can be read immidiatly as lists and since the output is a simple matter of expanding macros--to a certain extend, anyway), but the analysis of the program would still need to be done. Maybe it could be done more elegantly in lisp than C++, but for me it would have the additional overhead of having to familiarise myself with a new language.

/mailund

Just to continue the kibbitzing, if I were extending my version with the array representation optimizations you mention, I think I'd use an open implementation with optional representation hints, as I mentioned above for the FIB function, and worry about analyses to produce the hints automatically (and statically ensure correctness) as the very last step, that might not turn out to be needed in practice. But I don't intend to do the work, and don't know anything about the algorithms you're concerned with, so this opinion is only a wild guess.

It would be very interesting if you could post a simple algorithm that is representative of the ones that require the full power of DPROG -- if any simple one exists.

It could be very interesting to extend Darius's system to support such an algorithm and do the pruning of old values. I've often heard of C++ being favoured over Lisp in these applications for performance reasons, and I'm very curious to measure the actual performance difference on a realistic example. I don't know off hand whether to put my money on a factor of two, ten, or a hundred. Enquiring minds want to know :-)

The affine gap cost alignment example in the tutorial is actually such an example. (I'm looking at somewhat more complicated algorithms but the main difference is more tables and more bookkeeping--the essential stuff is found in the example).

There is three matrices that refer to each other, but in such a way that it is never necessary to store more than two rows of each to calculate all the cells, thus reducing the memory usage from O(N^2) to O(N) where N is the length of the sequences (which can easily be in the tens or hundreds of thousands). (In this example you can save a little more memory than needed to store the three times two rows, but you will still have linear memory usage so I don't bother with it).

In the general case you will have more matrices and you will need to store k rows for each matrix (where k may vary with the matrix as the rows in one matrix depend on the needs of other matrices).

It is not really that hard to calculate the number of rows I need to store for each matrix, it is just cumbersome to have to do it every time I change the algorithm. My hope is, that through DPROG it will be easier to prototype algorithms, and that the performance of the generated code is not too far from hand-optimized code. (At this point I have no idea about the performance of the generated code--I don't want to consider performance before I have the functionality).

As for relative performance between the generated C++ and any lisp implementation, I have no idea. I wouldn't be supprised if they were in the same ballpark if the lisp code is compiled. But, as I said, I have no idea.

/mailund

Python's re module uses this system, for example.

But I need to make sure that I never, ever, have a cache miss. If I have a cache miss I need to recalculate the missing value. Potentially this could change the complexity of the algorithm from, say, squared or cubed to exponential.

It is not enough to know the bound on cells I need to store, I also need to known which cells to store. If I know that, there is no need for the memorization approach; I can use dynamic programming.

This results in very simple code --- hardly more than the memoizer implementations already posted --- and often speeds your program up enough that you don't need more sophisticated garbage collection algorithms, like LRU or algorithmic analysis to determine what your program will need in the future.

I don't really need any fancy GC algorithm since I already know which values will be used and which will not from the program analysis. So I don't really use garbage collection, I just index into tables modulo some known bound.

/mailund