Gibbs sampling in lisp compared with C

I was reading two posts on Gibbs sampling: one by Daren Wilkinson in which he compares the performance of different programming languages on the same problem, and the other by Douglas Bates in which he performs the same task in Julia with impressive speed gains. For few weeks now, I’ve wanted to experiment with gsll (the GNU Scientific Library for lisp), specifically with their random number generators (by the way, here is a nice post by Jorge Tavares on the subject of generating random numbers using gsll). I thought doing the same benchmark with lisp might be good idea.

The reference C code

I will lift the C code from Daren Wilkinson:

 #include <stdio.h>
 #include <math.h>
 #include <stdlib.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>

 void main()
 {
   int N=50000;
   int thin=1000;
   int i,j;
   gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
   double x=0;
   double y=0;
   printf("Iter x y\n");
   for (i=0;i<N;i++) {
     for (j=0;j<thin;j++) {
       x=gsl_ran_gamma(r,3.0,1.0/(y*y+4));
       y=1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
     }
     printf("%d %f %f\n",i,x,y);
   }
 }    

I compiled and ran the code as suggested.

 gcc -O4 -lm -lgsl -lgslcblas gibbs.c -o gibbs
 time ./gibbs > datac.tab

As a matter of fact I ran ten times. The average came about 11.5 seconds on my Debian box at work.

The lisp code

You will need the gsll library which is on quicklisp. Get it and load it up.

 (ql:quickload :gsll)

Here is Wilkinson’s code translated (with a slight modification) into Common Lisp:

 (time 
    (with-open-file (res "datalisp.tab" 
                         :direction :output 
                         :if-exists :supersede)
      (let ((N 50000)
            (thin 1000)
            (rng (gsll:make-random-number-generator gsll:+mt19937+))
            (x 1.0d0)
            (y 0.0d0))
        (dotimes (i N) 
           (loop repeat thin do
               (setf x (/ (+ 1.0d0 (gsll:sample rng :gamma 
                                                     :a 3.0d0 
                                                     :b (/ (+ 4.0d0 (* y y)))))))
               (setf y (+ x (gsll:sample rng :gaussian 
                                             :sigma (* (sqrt 5.0d-1) x)))))
           (format res "~5D ~8,6F ~8,6F" i x y)))))

I ran the code 5 times on sbcl on the same Debian box at work. The average was about 57.5 seconds. This makes the common lisp code better than python but worse than PyPy in Wilkinson’s scale.