[PYTHON MATRIX-SIG] Alpha Matrix Object in C available - Guinea Pigs wanted

[James Hugunin]

I finally have an alpha version of a working, reasonably efficient matrix  
object implemented in C.  It compiles sucessfully on a Sparc10 running  
SUNOS, and on a P5(and P6) running NeXTStep.

I think that this object should satisfy many of the features requested by  
members of this group.

1) Extremely fast basic arithmetic functions

To me, this is sine qua non of a useful matrix object.  If I can't do  
arithmetic at close to the speed of hand coded C, then I have no use for the  
object.  As a sanity check that things are reasonably efficient, I  
performed a very rough comparision of the basic speeds of matlab, octave,  
python, and C for vector arithmetic.

The operation was to multiply a 10000 length vector of doubles with itself  
1000 times (10 M floating point multiplies).  The tests were all run on an  
unloaded Sun Sparc 10.

tool						speed(in MFlops)
Array.py in python				0.002
for loop in python				0.03
octave						0.2
matlab						1.6
matrixobject in python w/-O 			2.1	
C w/-O						2.4
C w/-O4						2.6

So, the python object is within 10% of the speed of the hand-coded C.  This  
is good enough for me.

2) Fancy arithmetic operations borrowing concepts from J.  Every operator  
can be given a rank, and can be used to perform direct, outer, and inner  
products as well as reductions and accumulations.  Much of the style of this  
part of the code is borrowed from Konrad's Array.py object.  Thanks Konrad  
for the introduction to J!

3) Very general multidimensional slicing operation.  This is based on Jim  
Fulton's proposal for multidimensional slicing, which those of you who have  
been following this list have seen.  This is a superset of most other  
slicing approaches that I've seen.  (Though I think that I need to go over  
Chris Chase's post on slices a little bit more carefully).

4) reshaping, general transposition, byteswaping, interface to/from strings, ...

What's obviously missing:

1) Complex numbers don't work right yet.  My basic problem is that I need a  
good complex number class in C to interface with, and I just haven't felt  
like writing this myself yet.

2) outer, inner, reduce, and accumulate style arithmetic functions could be  
optimized by a factor of 2-4 for simple cases.

3) inner product specfication is not completely clear to me when applied to  
operators with non-zero rank.

4) Documentation and comments in the code are minimal.

5) I'm sure that there must be a couple of memory leaks left hanging around.

What remains to be argued about (in my opinion at least)

1) Matrix-style multiplication as an option.  I really hate to leave the  
linear-algebra folks out in the cold.  On the other hand, I don't like the  
idea of setting a flag to determine how the "*" operator will act on any  
given object.  I'm open to suggestions.

2) Return by-value vs. by-reference.  I finally made the decision that  
every matrix slicing function returns by-value, and that indexing by a  
single value returns by-reference.  The by-value decision was based on a  
couple of bad experiences with the complexity of doing by-reference returns  
"right", and my observation that even in code that relied heavily on slicing  
a matrix I couldn't identify a performance difference that was ever greater  
than 20%.  Having a single index return by-reference is necessary to allow  
typical python style "multi-dimensional" indexing (ie. a[0][0][0]) to remain  
efficient (this great hack is Jim Fulton's idea).

3) Type-coercion.  At the moment, if I try to add a matrix of ints to a  
matrix of floats, I'll get an exception.  I'm not at all sure that it  
wouldn't be better to just coerce the matrix of ints to floats and then do  
the add.

4) Default type of a matrix.  At the moment, if I say "Matrix(1,2,3)", I'll  
get back a matrix of doubles.  This should probably be setable with a  
matrix module function, so that if I'm working a lot with complex numbers, I  
can make this return a complex matrix by default instead, or whatever other  
type I like to work in.

5) Rank-system.  I've really grown to like the J-style rank operators.  One  
thing that they provide is that I can say
Matrix(1,2,3)+Matrix([1,2,3],[11,12,13]) -> Matrix([2,4,6], [12,14,16]).
There are those who might believe that it would be safer for this to return  
an error that the two matrices aren't the same size.

6) There are a couple of possible patches to "ceval.c" to modify python to  
deal a little bit better with an object that is a sequence type, but that  
doesn't implement the sequence copy method when multiplied by a scalar.   
These would general-purpose changes and should be completely compatible with  
all existing python code.  At the moment, I've left this patch out.

7) Many other things that I'm sure I'm missing, which leads me to...

I'd really like to have people play around with this object a bit and offer  
me as harsh feedback as you wish.  This code is all based on Jim Fulton's  
implementation of a matrix object for using to interface to FORTRAN  
libraries,  Unfortunately, his employer has a small restriction on the  
distribution of this code (until it becomes officially released).

    THIS SOFTWARE HAS *NOT* BEEN APPROVED FOR RELEASE TO THE PUBLIC.

    THIS SOFTWARE MAY ONLY BE PROVIDED TO INDIVIDUALS WHO ARE
    CO-DEVELOPERS, REVIEWERS, OR TESTERS OF THE SOFTWARE, OR TO
    PERSONNEL OF THE U.S. GEOLOGICAL SURVEY.

So, anybody who wants to play around with this object, send me a message  
(at hugunin@mit.edu) saying that you agree to test and review it, I'll then  
point you to the appropriate FTP site.

If I get some positive feedback from people saying that they like the  
overall design of the object, then I'll go in and write up documentation and  
comment the code better.  'til then I'm not sure that I wouldn't just be  
wasting my time.

-Jim

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

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