[Numpy-discussion] What's wrong with matrices?

[Travis Oliphant]

Bill Baxter wrote:
> On 7/7/06, *Ed Schofield* <schofield at ftw.at <mailto:schofield at ftw.at>> 
> wrote:
>
>     Okay ... <Ed rolls up his sleeves> ... let's make this the thread ;)
>     I'd like to know why you, Sven, and anyone else on the list have gone
>     back to using arrays after trying matrices.  What was inconvenient
>     about
>     them?  I'd like a nice juicy list.  The whole purpose of the matrix
>     class is to simplify 2d linear algebra.  Where is it failing?
>
>
> Okay, here are a few that come to mind.
> 1) Functions that take a matrix but return an array.  Maybe these are 
> all fixed now.  But they better be fixed not just in numpy but in 
> scipy too.  To me this implies there needs to be some standard idiom 
> for how to write a generic array-protocol-using function so that you 
> don't have to think about it. 

A lot of these are fixed.  The mechanism for handling this is in-place: 
either using asanyarray in the function or (more generally) using a 
decorator that wraps the arguments with asarray and returns the output 
with __array_wrap__.

But, we need people to help with fleshing it out.   

>
> 2) At the time I was using matrix, scalar * matrix was broken.  Fixed 
> now, but that kind of thing just shouldn't happen.  There should be a 
> tests for basic operations like that if there aren't already.

We need people to write and implement the tests.   It's one way 
everybody can contribute.   I do use matrices occasionally (not never as 
has been implied).  But, I do more coding than linear algebra 
(particularly with images), therefore my need for matrix math is 
reduced.   Nonetheless, I've been very willing to make well-defined 
changes that are needed. 

Most problems cannot be found out without people who use and test the 
code. 
>
> 3) mat() doesn't make sense as a shortcut for matrix construction.  It 
> only saves 3 letters over typing matrix(), and asmatrix is generally 
> more useful.  So mat() should be a synonym for asmatrix()
I'd be willing to make that change, but it will break some people's 
SciPy code.
>
> 4) eye,empty,rand,ones,zeros,arange and anything else that builds an 
> array from scratch or from a python list should have a matrix equivalent
Would

from numpy.defmatrix import ones, zeros, ...

work?
>
> 5) I've got squeezes like crazy all over my matrix-using code.  Maybe 
> this was a bug in 0.9.5 or so that's been fixed?  I do seem to recall 
> some problem with indexing or c_ or something that was causing 
> matrices to grow extra levels of length 1 axes.  Again, like the 
> scalar*matrix bug, things like that shouldn't happen.
Sure, but it's going to happen in a beta release...  That's why we need 
testers.   As I recall, most bugs with matrices have been fixed fairly 
quickly as soon as they are reported.

>
> 6) No good way to do elementwise operations?  Sometimes you just want 
> to do an elementwise mult or divide or exponentiation.  I guess you're 
> supposed to do  Z = asmatrix(X.A * Y.A).  Yah, right.
This is a problem with a dearth of infix operators.   In fact, if we had 
a good way to write matrix multiplication as an infix operator, perhaps 
there wouldn't be any need for matrices. 

I'm really not sure how to fix the problem (the .M attribute of arrays 
was an attempt to make it easier):

(X.A * Y.A).M

But, there is always

multiply(X,Y)

>
> 7) Finally, once all that is fixed, I find the slavish adherance to 
> ndim=2 to be too restrictive. 
>     a) Sometimes it's useful to have a stack of matrices.  Pretty 
> often, in fact, for me.  I guess I could use a python list of matrices 
> or an object array of matrix or something, but I also think there are 
> times when it's useful to treat different pairs of axes as the 
> 'matrix' part.   So I'd like matrices to be able to have ndim>2.  
I suppose this restriction could be lifted. 

>     b) On the other end, I think ndim<2 is useful sometimes too.  Take 
> a function like mean(), for example.  With no arguments  the return 
> value is a 1x1 matrix (as opposed to a scalar).
Have you checked lately.  It's a scalar now...  This has been fixed.

> Or take indexing.  It seems odd to me that where() reurns a tuple of 
> shape==(1,N) objects instead of just (N,) . 
The way to fix some of these is to return arrays for indexing instead of 
allowing matrices.    But, matrices that are less than 2-d just don't 
make sense. 
> Maybe I can get over that though, as long as it works for indexing 
> (which it seems it does).   But I think the scalar return case is a 
> real issue.  Here's another: sum().  For an array you can do 
> sum(sum(a)) and get a scalar if a is 2-d, but for matrix sum(sum(m)) 
> is the same as sum(m).   And along these lines, m[newaxis] just 
> silently doesn't do anything.  That doesn't seem right.

These are just semantic questions.   It's no surprise that sum(sum(m)) 
returns the same as sum(m) for a matrix because summing over the same 
axis won't change the result.   You have to sum over both axes in a matrix.

Thanks for the feedback.

-Travis