A Bet Over Beer (or how to write a faster root() function in IML in SAS 9.4)

I haven't been here in a while -- previous job did not allow blogging.  

Over beers the other night with some friends and alumni from SAS, someone mentioned IML is now multithreaded. Someone else pipes up with, "yeah but the Cholesky root() function is not and that's all I really use it for."

To which I respond, "A -- I didn't know that about multithreading in  SAS 9.4.  Very cool.  B -- I bet I can write a faster root() function for you."  What hubris!  I blame the beer.

So having thrown down the gauntlet, I went back to my hotel room, dug out old posts (here), and got to work.

Here's what I came up with:

start my_chol(a,b);

       n = nrow(a);

       out = a;

       x = floor(n/b-1)*b;

       i=1;

       do while(i<x);

             out[i:(i+b-1),i:(i+b-1)] = root(out[i:(i+b-1),i:(i+b-1)])`;

             out[(i+b):n,i:(i+b-1)] = out[(i+b):n,i:(i+b-1)] * inv((out[i:(i+b-1),i:(i+b-1)])`);

           out[(i+b):n,(i+b):n] = out[(i+b):n,(i+b):n] - out[(i+b):n,i:(i+b-1)]*(out[(i+b):n,i:(i+b-1)])`;

             i = i + b;

       end;

       out[(x+1):n,(x+1):n] =root(out[(x+1):n,(x+1):n])`;

       do i=1 to n-1;

             out[i,(i+1):n] = J(1,n-i,0);

       end;

       return(out);

finish my_chol;

This takes 2 parameters.

1. A – The matrix to be decomposed.
2. B – A block size.  This trades off the multithreading in matrix multiplication with single threaded root() and inv() functions. You need to play with this to see what works best on your system. 100-150 seems to be the sweet spot for my laptop.

Running this on a 5000x5000 matrix (B=100), I get:

1. IML root() – 21.08 seconds
2. my_chol() – 13.98 seconds

I’ll collect my free beer later this week J.