[Numpy-discussion] Matrix Class [was numpy release]

[Charles R Harris]

On Thu, Apr 24, 2008 at 5:45 PM, Timothy Hochberg <tim.hochberg at ieee.org>
wrote:

>
>
> Chris.Barker wrote:
>
>> Alan G Isaac wrote:
>> > the cost of complexity should be justified by a gain in functionality.
>>
>> I don't think functionality is the right word here. the Matrix class(es)
>> is all about clean, convenient API, i.e. style, not functionality -- we
>> have all the functionality already, indeed we have it with plain old
>> arrays, so I think that's really beside the point.
>
>
> Not entirely, there's no good way do deal with arrays of matrices at
> present. This could be fixed by tweaking dot, but it could also be part of a
> reform of the matrix class.
>
>  [CHOP]
>
>>
>> Timothy Hochberg wrote:
>> >    1. The matrices and arrays should become more alike if possible
>>
>> I'm not sure I agree -- the more alike they are, the less point there is
>> to having them.
>
>
> That's the best possible outcome. If some solution can be reached that
> naturally supports enough matrix operations on array, without significantly
> complexifying array,  to satisfy the matrix users then that would be great.
> Less stuff to maintain, less stuff to learn, etc, etc.
>
> With that in mind, what is minimum amount of stuff that matrix should
> support:
>
>    1. Indexing along a single axis should produce a row or column vector
>    as appropriate.
>    2. '*' should mean matrix multiplication.
>    3. '**' should mean matrix exponentiation. I suspect that this is less
>    crucial than 2, but no more difficult.
>
> There's some other stuff that's less critical IMO (.H and .I for example)
> but feel free to yell at me if you think I'm mistaken.
>
> There's some other properties that a fix should have as well, in my opinion
> and in some cases in others opinions as well.
>
>    1. A single index into an array should yield a 1D object just as
>    indexing an array does. This does not have to inconsistent with #1 above;
>    Chris Barker proposed one solution. I'm not sold on the details of that
>    solution, but I approve of the direction that it points to. [In general
>    A[i][j] should equal A[i,j]. I know that fancy indexing breaks this; that
>    was a mistake IMO, but that's a discussion for another day].
>    2. It should be possible to embed both vectors and matrices naturally
>    into arrays. This would allow natural and efficient implementation of, for
>    example, rotating a 1000 3D vectors. One could spell that R * V, where R is
>    a 2x2 matrix and V is a 1000x3 array, where the last axis is understood to
>    be a vector.
>
>
R is 3x3? This works already, whats tough is an array of rotation matrices
applies "element wise" to an array of vectors.


>    1.
>    2. I'm pretty certain there's a third thing I wanted to mention but it
>    escapes me. It'll resurface at the most inopportune time....
>
> Let's play with Chris Barker's RowVector and ColVector proposal for a
> moment. Let's suppose that we have four four classes: Scalar, RowVector,
> ColVector and Matrix. They're all pretty much the same as today's array
> class except that:
>
>    1. They are displayed a little differently so that you can tell them
>    apart.
>    2. __mul__ and __pow__ are treated differently
>
> Let's consider __mul__: when a RowVector multiplied with ColVector, the dot
> product is computed. If, for example, the arrays are both 2D, the they are
> treated as arrays of vectors and the dot product of each pair of vectors is
> computed in turn: I think broadcasting should work correctly, but I haven't
> thought that all the way through yet. Ignoring broadcasting for a moment,
> the rules are:
>
>    1. (n)D Scalar * (n)D Scalar => (n)D Scalar (elementwise)
>    2. (n)D RowVector * (n)D ColVector => (n-1)D Scalar (dot product)
>    3. (n+1)D Matrix * (n)D ColVector => (n)D ColVector (matrix-vector
>    product)
>    4. (n)D Matrix * n(D) Matrix => (n)D Matrix (matrix) product
>
> And ColVector * RowVector is a matrix? Or illegal. The latter, I suppose,
if * is contraction on indices.


> Other combinations would be an error. In principal you could do dyadic
> products, but I suspect we'd be better off using a separate function for
> that since most of the times that would just indicate a mistake. Note that
> in this example Scalar is playing the role of the present day array, which
> I've assumed has magically vanished from the scene somehow.
>
> This looks like it gets of most the way towards where we want to be. There
> are some issues though. For example, all of the sudden you want to different
> transpose operators; one that transposes matrices and flips colVectors to
> RowVectors and leaves Scalars alone, and another that manipulates the rest
> of the array structure. There is probably other stuff like that too, it
> doesn't look insurmountable to me right now, however.
>
> Still, I'm not sold on the whole RowVector/ColVector/Matrix approach. I
> have a gut feeling that we'd be better off with a single array class that
> was somehow tagged to indicate it's type. Also, I keep hoping to come up
> with an elegant generalization to this that turns arrays into quasi-tensor
> objects where all the matrix behavior falls out naturally. Sadly, attempts
> in this direction keep collapsing under there own weight but I'm still
> thinking about it.
>

I think a start could be made by adding a signature to arrays, indicating
whether the corresponding indices are up or down. Up is column, down is row.
The in tensor parlanc, matrices have signature (u,d). The matrix product
sums an up with a down, contraction, so you get the usual rules. Note that
row*col == col*row == inner product in this scheme, which may seem a bit
odd, but transpose(row) := col as usual. This doesn't solve the problem of
containers of matrices, however. I think that matrices could be treated as
objects with special, efficient storage in arrays, and that way the
element-wise operators will do what you need using broadcasting. Essentially
the last two/one dimensions are reserved. This is a bit like matlab cell
arrays, except the arrays must be homogeneous in the matrix type. Note that
I'm making a distinction between 1xn matrices and row vectors. The former
has signature (u,d), the latter (u,).

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080424/88191ac0/attachment.html>