[PYTHON MATRIX-SIG] The Matrix Object Proposal (very long)
James Hugunin
jjh@ling-ling.lcs.mit.edu
18 Aug 1995 17:16:55 GMT
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
=================