[Numpy-discussion] Re: array vs. matrix for matlab-compatible stuff

[Schofield, E. <Ed....@ft...>]

(Apologies if this is a duplicate.  My mail isn't getting through ...)

Travis Oliphant wrote:
>> Sven Schreiber wrote:
>>
>> I see, thanks for the quick answer. So wouldn't it be a good idea to
>> have all the specifically
>> matrix-related stuff (afaics, things in numpy/lib/twodim_base.py)
>> return matrices?
>>
>> It seems my question (or misunderstanding) has a broader scope:
>> Coming from matrix languages, I'm
>> glad about short notations like A.I or A*B representing standard
>> matrix operations. Much easier than
>> linalg.inverse(A) or matrixmultiply(A,B). However, many matrix
>> functions (decompositions etc.) seem
>> to return arrays instead of matrices, even if you feed them matrices
>> (is this assumption correct?).
>> So you have to use mat(returned-result) again and again to be able to
>> do return-result.I afterwards,
>> which seems clumsy. So what's the best strategy here?
>
> Yes, some of them do still return arrays.  Matrices are longer lived
> in NumPy then the were in Numeric, for sure, but many functions still
> aren't friendly to matrices and convert all inputs to arrays before
> operation.   Originally, I had the asarray(...) function not convert
> matrices by default, but this is too surprising because matrices
> change the '*' and '**' operators which could make your function not
> work.
> We should convert all the functions that don't handle matrices so they
> will.  I'd like to see matrices survive longer than they do.   There
> are some functions that try to do that

I agree, and I'm willing to help with this.  We should also think about
handling SciPy's spmatrix objects properly.  In this case we don't want
NumPy to have any dependency on a particular format for SciPy's spmatrix
objects -- but we could design and publish a simple interface for sparse
matrix objects (from SciPy or anywhere else) to conform to for basic
NumPy compatibility.

How should we go about this?  Let's take linalg.solve_linear_equations()
as an example.  We could use isinstance(x, matrix) for dense matrix
objects, because they're part of NumPy, then use asarray(x), process =
them
as arrays, and wrap them up with asmatrix(x) at the end.  But with =
sparse
matrices this wouldn't work.  What if we specify instead that objects
passed to functions like solve_linear_equations() need to provide two
members:

    x.toarray(), allowing x to represent itself as an array
    x.fromarray(a), creating another object of the same type as x out of
                    the given array a

We could choose other names, of course, maybe asarray() or whatever =
else.
But this would simplify the code for all these functions while allowing
us to support other array-like objects without needing to know any more
about them.

I'm also +1 on eye() returning a matrix :)

-- Ed