[Numpy-discussion] doing zillions of 3x3 determinants fast

[Daniel Lenski]

Hi all,
I need to take the determinants of a large number of 3x3 matrices, in 
order to determine for each of N points, in which of M tetrahedral cells 
they lie.  I arrange the matrices in an ndarray of shape (N,M,5,3,3).

As far as I can tell, Numpy doesn't have a function to do determinants 
over a specific set of axes... so I can't say:
  dets = np.linalg.det(a, axis1=-1, axis2=-1)

Using an explicit Python loop is very slow... doing the determinants of 
100,000 random 3x3 matrices in a loop and inserting the results into a 
preallocated results array takes about 5.1 seconds.

So instead I coded up my own routine to do the determinants over the last 
2 axes, based on the naive textbook formula for a 3x3 determinant:
  def det3(ar):
    a=ar[...,0,0]; b=ar[...,0,1]; c=ar[...,0,2]
    d=ar[...,1,0]; e=ar[...,1,1]; f=ar[...,1,2]
    g=ar[...,2,0]; h=ar[...,2,1]; i=ar[...,2,2]
    return (a*e*i+b*f*g+c*d*h)-(g*e*c+h*f*a+i*d*b)

This is *much much* faster... it does the same 100,000 determinants in 
0.07 seconds.  And the results differ from linalg.det's results by less 
than 1 part in 10^15 on average (w/ a std dev of about 1 part in 10^13).

Does anyone know of any well-written routines to do lots and lots of 
determinants of a fixed size?  My routine has a couple problems I can see:
  * it's fast enough for 100,000 determinants, but it bogs due to 
    all the temporary arrays when I try to do 1,000,000 determinants
    (=72 MB array)
  * I've heard the numerical stability of the naive determinant 
    algorithms is fairly poor, so I'm reluctant to use this function on 
    real data.

Any advice will be appreciated!

Dan