[PYTHON MATRIX-SIG] The Matrix Object Proposal (very long)

[James Hugunin]

I apologize for the length of this posting, but you were warned in the  
subject header.

There seems to be a fair amount of interest in the python community  
concerning the addition of numeric operations to python.  My own desire is  
to have as large a library of matrix based functions available as possible  
(linear algebra, eigenfunctions, signal processing, statistics, etc.).  In  
order to ensure that all of these libraries interoperate, there needs to  
be agreement on a basic matrix object that can be used to represent arrays  
of numbers.  The following is a proposal for such an object.  It can be  
viewed as an extension of the current array object to higher dimensions  
and numerical operations.

The wheels have already been set in motion to create a Matrix SIG  
mailing-list for further discussion of this proposal, as well as other  
python numerics related issues.  This should come on-line sometime next  
week and serious discussion of this proposal should probably be deferred  
until it can be done on that list.  Summaries of that discussion will be  
posted to the main list as consensus develops.

Below I specify the behavior of a Matrix object within python.  The C (or  
FORTRAN) API to such an object is obviously of significant importance as  
this object is primarily intended as an interface to existing numerical  
libraries.  I'm currently working on that proposal (at least the C part)  
as well, but I think everyone will agree that this post is long enough as  
it is.

Note: Parts of this proposal are stolen from the python documentation for  
the array object, and from discussions on the net.  I'd like to  
particularly thank Jim Fulton for letting me review (and steal ideas from)  
his current Matrix object and for providing feedback on the first draft of  
this proposal (though he certainly doesn't agree with everything in it,  
even now).


Matrix Object

This defines a new object type which can efficiently represent a  
multidimensional array of basic values (char, unsigned byte, signed byte,  
int, short, long, float, double, complex float and complex double).   
Matrices are defined as both a sequence type and a mapping type, and  
usually (unless it's a matrix of char's) as a number type.  It should be  
noted that allowing a type to be both a sequence and a number and to  
behave properly for addition and multiplication, will require some very  
small patches to the ceval.c code.

All of the examples use the following matrices:
Assume A = [1,2,3], B = [11,12,13], C = [[1,4,9],[16,25,36]],
M=[[ 1, 2, 3, 4, 5],[11,12,13,14,15],[21,22,23,24,25],[31,32,33,34,35]]


**Numeric operations**

"+", "-", "*", "/", "^"(as power), abs, unary"-", and unary"+" are defined  
as the corresponding element-wise arithmetic operations.  Because these  
assume element-wise correspondence between the two matrices, the operandi  
must be of the same dimensions (however, note the dimensions-coercion  
section below).

A+B = [12,14,16]
A*B = [11,24,39]
C*C = [[2,8,18],[32,50,72]]
-A = [-1,-2,-3]

Possible suggestions for the other numeric operations:
"%" represents matrix multiplication (or vector dot product).
"~" represents transposition.
Other suggestions?


**Sequence operations**

Concatenation and multiplication sequence operations are not defined, but  
instead the arithmetic operations are invoked for the "+" and "*"  
operators.  Indexing by both single indices and by slices is supported.   
All returned matrices are returned by reference rather than by value.
len(M) is defined to return the size of the first dimension in in matrix.


examples:
len(M) -> 4
 
D = M[0]
D[0] = 66
print M[0]
	[66,2,3,4,5]

The semantics of returning a slice by reference is consistent with the  
mapping operations given below, however, this is inconsistent with the  
list objects semantics.  This is open to debate.

D = M[0:2]
D[0][0] = 66
print M[0]
	[66,2,3,4,5]


**Mapping operations**

Recent newsgroup discussions have offered some convincing arguments for  
treating a matrix as a mapping type where sequences of ints are used as  
the keys.  The sequence must be of length less than or equal to the numer  
of dimensions in the matrix.  Each element of the sequence is either an  
integer or a sequence.  If it is an integer, it returns the corresponding  
element for that dimension, if it is a sequence then it returns all of the  
elements given in the sequence.

ie. A[0] -> 1, A[((1,2))] -> [2,3]
M[0] -> [1,2,3,4,5]
M[(0,3)] -> 4
M[((0,2),0)] -> [1,21]
M[(range(1,3),range(2,4))] -> [[13,14],[23,24]]

In order to simplify the syntax of these references, it would be valuable  
to change the syntax of python so that M[0,3] is equivalent to M[(0,3)].   
This is a small change to the current indexing semantics which will give a  
reasonable meaning to a currently illegal syntactic item.

When used as a setvalue, the corresponding elements of the given array  
will be set to the left hand side.  All of these operations that return a  
matrix will return it by-reference, meaning that the elements of the new  
array will be references to the old one and changes to either array will  
effect the other.  


**Instantiation**

The following function will be used to instantiate a matrix:
Matrix(datasource(optional), d1, d2, ..., typecode='d')

The typecodes are as follows:
c - char (non-numeric)
1 - unsigned char
b - signed char
h - short
i - int
l - long
f - float
d - double
F - complex float
D - complex double

The optional datasource can be a file object or a string object and will  
fill the matrix with the raw data contained therein.  In this case the  
dimensions are required.  The first dimension can be input as -1, in which  
case the first dimension will be set to be as large as necessary to hold  
the entire contents of the file or string.

The optional datasource can also be a sequence, in which case the  
dimensions of the array do not need to be passed in and are instead  
determined by the hierarchical structure of the sequence.

examples:
Matrix([range(3),(4,5,6)], 'i') -> [[0,1,2],[4,5,6]]
Matrix(2,3) -> [[0,0,0],[0,0,0]]
Matrix('abcdef', 2,3, 'u') -> [[97,98,99],[100,101,102]]
Matrix('abcdef', -1,2, 'c') -> ['ab', 'cd', 'ef']


**Coercion**

Dimension coercion:The following is a proposal to make the matrix object  
behave as closely as possible to matrices in matlab (to make it easy to  
convert matlab code to python), and to generalize this behavior to  
arbitrarily high dimensional objects.

If a matrix of fewer dimensions than that required is provided, then the  
matrix can be converted to more dimensions by making copies of itself.   
For purposes of dimension coercion, an int or a float is considered a 0-d  
matrix.

ie. A+1 -> [2,3,4], A+C -> [[2,6,12],[17,27,39]]
	B^2 -> [121,144,169] vs. B^A -> [11, 144, 2197]

Simliar coercion should be performed for most functions acting on  
matrices.  Helper functions will be provided to make this easy, however  
there is no way to enforce this.  The following is a proposed standard to  
allow functions to be applied to higher-dimensional objects in a uniform  
way.

If a function is defined to operate on scalars (say sqrt()), then when  
applied to a matrix it should apply itself to every element in the matrix.

ie. sqrt(C) -> [[1,2,3],[4,5,6]]

This is generalized to higher dimensions so that if a function is provided  
with a matrix of more dimensions that it requires, that function will  
apply itself to each appropriately sized unit in the matrix.

ie. sum(A) -> 6 whereas sum(C) -> [14,77]


Type coercion: No type coercion will be performed.  However, it is  
possible to have a function defined to work on multiple types (operator  
overloading).  So sqrt might be defined on matrices of both floats and  
doubles in which case it will work properly when given either matrix as  
input, but will return a type error if applied to a matrix of ints.


**Basic Data Items and Methods (non-numeric)**

typecode - The typecode character used to create the matrix

itemsize - The length in bytes of one matrix item in the internal  
representation

dimensions - A sequence containing the size of the matrix along each of  
its dimensions

copy - Returns a copy of the matrix.  This is a good way to take an array  
returned by reference and turn it into one returned by value.

transpose - The transpose of the matrix (not needed if "~" is used for  
transpose)

byteswap - Useful for reading data from a file written on a machine with a  
different byte order

typecast(t) - Returns a new matrix of the same dimensions where each item  
is cast (using C typecasting rules) from the matrix's type to the new type  
t.

tofile(f) - Write all items (as machine values) to the file object f.

tolist() - Convert the matrix to an ordinary (possibly nested) list with  
the same items.

tostring() - Convert the array to an array of machine values and return  
the
string representation (the same sequence of bytes that would
be written to a file by the tofile() method.)

Note: the append, fromlist, etc. methods are not present because it is  
assumed that a matrix object is of constant size (there are some strong  
efficiency reasons for this assumption but this is likely to be another  
source of arguments).

**Basic functions taking matrix arguments (non-numeric)**

Since it isn't possible to append to a matrix (or efficient if it's  
decided it should be possible) I thought that some sort of function for  
producing a new matrix from a list of matrices would be useful.  So I'm  
proposing the following function which is very similar in spirit to the  
string.join operation.

concat(l) - l must be a list of matrices, and each matrix in l must have  
the same number of dimensions, as well as the same size for every  
dimension except for the first.  This will return a new matrix M whose  
first dimension is the sum of the first dimensions of all the matrices in  
the list, and whose data is filled from the data in all of these matrices.

ie. concat([A,B]) -> [1,2,3,11,12,13]

This is different from Matrix([A,B]) -> [[1,2,3],[11,12,13]]


**Basic functions taking matrix arguments (numeric)**

It is hoped that there will be a large variety of special purpose modules  
to do fancy numeric operations on matrices.  There should also be a basic  
set of operations provided in the Matrix module (the module with the  
instantiation function).  Because type coercion is not going to be  
performed, it is important to decide both on the list of functions and on  
the basic types they should support.  The following is a very preliminary  
list of the functions I consider most fundamental.

rng(start, stop, step(optional), typecode='d') where all points in the  
range must be strictly less than the stop point, according to the  
semantics of range().  This function should possibly be added to the  
instantiation options.

ie. rng(0,0.4) = [0,0.1,0.2,0.3]; rng(0,0.1,0.41) = [0,0.1,0.2,0.3,0.4]


scalar operations (defined on floats, doubles, and all complex):
everything in mathmodule.c
round

scalar operations (defined on all numeric types):
sign

vector to scalar operations (defined on all numeric types):
max, min - Should these return the index of the max or min?
sum, prod

vector to vector operations (defined on all numeric types):
cumsum, cumprod

-----
Jim Hugunin - hugunin@mit.edu
Laboratory for Computer Science
Massachusetts Institute of Technology

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================