[Numpy-discussion] Stacking matrices

[Tim Hochberg]

The recent message by Ferenc.Pintye (how does one pronounce that BTW) 
reminded me of something I've been meaning to discuss: I think we can do 
a better job dealing with stacked matrices. By stacked matrices I mean 3 
(or more) dimensional arrays where the last two dimensions are 
considered to be matrices. Let us suppose that one wants to find the 
eigenvalues of 200,000 3x3 matrices as Ferenc does. At present you have 
to call linalg.eigvals 200,000 time; the overhead is going to kill you. 
I've run into this problem before myself in the past (and kludged my way 
around it), so the problem is not limited to Ferenc.

I see no reason that this could not be fixed for the linalg functions. 
The details would vary from function to function, but for eigvals, for 
example, when passed an NxMxM array, it would return an NxM array of  
eigenvalues. One can actually get pretty far just modifying the python 
functions in linalg.py, or at least one could in numarray and I'm 
guessing that numpy is similar. That is because there is a fair amount 
of overhead in setting up the function calls and this can all be done 
just once. I have code lying around that does this for solve, inverse 
and determinant that works with numarray that would probably not be hard 
to adapt to numpy if we were interested in making this change of 
interface to linalg.

Another place that this same effect has bitten me in the past is with 
dot. Dot does work on higher dimensional arrays, but it attempts to do a 
tensor product of sorts. While this sounds nifty, in practice I've found 
it less than useful. However, because dot computes both scalar and 
matrix products it's not clear that a set of rules that was more or less 
backwards compatible and that extended to higher dimensions in a more 
useful way than at present could be crafted. Probably something with 
different behavior, perhaps that always performed matrix products on the 
last two axes, could be added with a name other than dot.  This avoids 
backwards compatibility issues as well.

None of this stuff is worth holding up the the release 1.0 over, but I 
though I would mention it while it was on my mind.

-tim