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
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
