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

Daniel Lenski dlenski@gmail....
Sun Aug 24 22:48:54 CDT 2008


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



More information about the Numpy-discussion mailing list