[PYTHON MATRIX-SIG] Release 0.20 of matrix object is available

[Hinsen Konrad]

   Major incompatible changes (I expect that these might elict some  
   interesting discussion):

OK, let's start ;-)

   1) Matrix indexing now returns by reference to the original data, not by  
   value.  As a side effect of this change, arbitrary sequences are not allowed  
   in multidimensional indices, but only single indices, slices, RubberIndex  
   and None.

I don't see any problem with this. Of course the documentation (to be written...)
should be clear so that everyone knows the implications of this, but since
Python variables in general hold references, it shouldn't be a surprise
for anyone.

   This change was motivated both by the possible speed increases (about 40%  
   for some code) and more importantly by issues of clarity of expression.  I  
   couldn't come up with any other way to make the following hold:

   m[0][1] is an efficient way to index a multi-dimensional array.
   m[0,1] == m[0][1]
   m[0,] == m[0]

I don't see why these are important features, but of course they don't
do any harm.

   2) No more ranks of operators, instead Yorick style pseudo indices are  
   used.  As a side effect of this, outer product is now a convenience function  
   rather than a method on ofuncs.

As another side effect, we have lost some important functionality.

   For functions of unbounded rank (currently the only kind of function  
   supported by my omathmodule) Yorick pseudo indices support a superset of  
   what was possible to express using ranks and outer products.  There is no  
   fundamental reason not to have both pseudo indices and ranks, I just think  
   that it's conceptually cleaner not to mix the two up.

The pseudo indices are not a replacement for ranks, but just a convenient
notation for some structural functions. I like this notation very much,
but we still need ranks. And as you said, there is no reason not to have both.

Apart from the problem of functions with finite rank (which we will
certainly need e.g. for matrix inversion), it is not true that pseudo
indiced support a superset of what is possible with ranks and the
associated structural functions. In fact, it is just the opposite, as
I will demonstrate.

First let me explain how the new pseudo indices fit into the J-style
array system. J doesn't have the concept of "indexing"; instead there
are a few structural functions that can be applied with different
ranks to construct all kinds of rearrangements of array elements. All
the (pseudo) indices we currently have can be mapped onto these
structural functions as follows:

- Explicit numerical indices, lists, slices and the rubber index
  are equivalent to the function take() (see my Array.py) with
  various rank.
- The pseudo-index "None" is equivalent to the J-function "itemize"
  that adds an axis of length one to its argument. I didn't implement
  this in Array.py, but it is a trivial addition.
- The contracting rubber index (* in Yorick, I haven't found it yet
  in Python)is equivalent to the function ravel() with appropriate
  rank.

As an example, the outer product mentioned above could be
(inefficiently) implemented in my Array.py as:

def outer(f, a, b):
   return f(reshape[0](a, shape(b)), b)

So every construction with (pseudo) indices can be emulated with three
structural functions of suitable rank. On the other hand, there are
constructions that cannot be expressed with pseudo indices, for
example all cases where a function rank is not a constant but a
variable. In fact these can be handled in Python -- as opposed to
Yorick -- by remembering that a[x,y,z] is really just a shorthand for
a[(x,y,z)] and that the index tuple can be constructed from an
expression, but that seems more like an implementation accident than a
feature to me.

More importantly, pseudo indices don't help in the least when it comes
to reductions and general inner products that require the specification
of an axis. Yorick does this by putting the function into the index
list -- a clever idea, but I'd prefer to keep a clear distinction
between functions and arguments. The current matrix module allows
a second parameter to reduce() that indicates the axis, which is
exactly how APL handles the problem.

But both these methods are restrictive. They are sufficient only
for reductions with functions of unbounded rank. I can give an
example of an actual application where this is not enough: I had
a sequence of n rotation matrices, stored in a rank-3 array D of
shape (n, 3, 3). Then the net rotation of this sequence could
be obtained by matrixMultiply.reduce(D). Since matrix
multiplication is a rank-2 function, I couldn't express this
with the current matrix module.


The bottom line: pseudo indices are convenient shorthands for
the most common structural functions, but they are no replacement
for the general function rank scheme we had in version 0.11.


And now for something completely different
------------------------------------------

1) Names

This may seem to be a minor issue, but we should tackle this before
we all get used to the current ones.

I strongly dislike the type-specific constructors (also used for
output). They should be IntegerMatrix, FloatMatrix, ComplexMatrix,
CharacterMatrix, and GeneralMatrix. Anyone is free to define shorter
names for efficient typing, if desired.

I even propose a more radical renaming. Many people associate "matrix"
with the 2D-matrices from linear algebra. So it would be better to
call our general objects "arrays", and leave the name "matrix" for
linear-algebra type objects that are restricted to rank 2 and use *
for matrix multiplication.  They could be implemented in Python based
on arrays.

The pseudo-index 'None' is very confusing. An index "2" picks item
number two from an axis, an index "All" picks all items from an
axis. Consequently "None" should pick no item from an axis, which is
course is a pointless operation. What it really does is create a new
axis, so it should be named "New".  "RubberIndex" is not confusing,
but a bit strange, so it's worth thinking about alternatives. How
about "Skip", in the sense "skip all following axes for explicit
consideration"? For the "*" in Yorick I propose "Contract".


2) Constructors

Currently there are two constructors, matrix() and Matrix(), with
slightly different behaviour: matrix() takes a single argument that is
a (possibly nested) Python sequence object, e.g. a list. Matrix()
takes an arbitrary number of arguments that are made into a list
before being passed to matrix(). So Matrix(1,2,3) is equivalent
to matrix([1,2,3]).

I find it confusing to have to constructors with almost the same name
and almost the same function. I propose to have only one, of whatever
name, which behaves like matrix(). The reason is that many matrix
functions accept nested lists or tuples instead of matrix arguments,
e.g. matrix([2,3]) + [5,4] works. So does matrix([2,3]) +
matrix([5,4]), which is equivalent. But Matrix([2,3]) + [5,4] and
Matrix([2,3]) + Matrix([5,4]) are not equivalent (well, they are in
this simple exaple due to broadcasting, but not in general). So
matrix() can be considered a conversion function from lists and tuples
to arrays, which is needed anyway.

I realize that the calling style of Matrix() has the advantage of
saving one pair of brackets, but I don't consider this important
enough to create the current confusion.


I guess that's enough for today ;-) I wish the lucky participants
of the workshop much fun and hope that they don't have to eat
spam for lunch ;-)

-------------------------------------------------------------------------------
Konrad Hinsen                     | E-Mail: hinsenk@ere.umontreal.ca
Departement de chimie             | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. Centre-Ville     | Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase experimentale)
-------------------------------------------------------------------------------

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

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