[Numpy-discussion] Stacking matrices
tim.hochberg at cox.net
Thu Jul 20 17:07:24 CDT 2006
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.
More information about the Numpy-discussion