[PYTHON MATRIX-SIG] Alpha Matrix Object in C available - Guinea Pigs wanted
James Hugunin
jjh@mama-bear.lcs.mit.edu
Mon, 30 Oct 95 12:28:37 -0500
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
=================