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

Anne Archibald peridot.faceted@gmail....
Sun Aug 24 23:53:33 CDT 2008

2008/8/25 Daniel Lenski <dlenski@gmail.com>:
> On Mon, 25 Aug 2008 03:48:54 +0000, Daniel Lenski wrote:
>>   * 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 managed to reduce the memory usage significantly by getting the
> number of temporary arrays down to exactly two:
> 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]
>    t=a.copy(); t*=e; t*=i; tot =t
>    t=b.copy(); t*=f; t*=g; tot+=t
>    t=c.copy(); t*=d; t*=h; tot+=t
>    t=g.copy(); t*=e; t*=c; tot-=t
>    t=h.copy(); t*=f; t*=a; tot-=t
>    t=i.copy(); t*=d; t*=b; tot-=t
>    return tot
> Now it runs very fast with 1,000,000 determinants to do (<10X the time
> required to do 100,000) but I'm still worried about the stability with
> real-world data.

As far as I know, that's the best you can do (though come to think of
it, a 3D determinant is a cross-product followed by a dot product, so
you might be able to cheat and use some built-in routines). It's a
current limitation of numpy that there is not much support for doing
many linear algebra problems. We do have perennial discussions about
improving support for array-of-matrices calculations, and in fact
currently in numpy SVN is code to allow C code doing a 3x3 determinant
to be easily made into a "generalized ufunc" that does ufunc-like
broadcasting and iteration over all but the last two axes.


More information about the Numpy-discussion mailing list