[PYTHON MATRIX-SIG] Matrix multiplication

[carlos fonseca]

Hello,

I may have missed any previous discussion on the implementation of
matrix multiplication in Numerical Python. Anyway, I understand that
the current implementation of .matrixMultiply() requires arrays to be
2-dimensional (has this changed?).

While trying to port some matlab code to Numerical Python and make it
more general (eliminate for loops, add functionality, etc), I realised
that it would be nice to be able to:

1) treat (e.g.) the last 2 dimensions of an array as matrices, i.e.,
  having arrays of matrices. For example,

  A[m-by-n-by-p] * B[p-by-q] = A[m-by-n-by-p] * B[1-p-by-q]
                             = C[m-by-n-by-q]

  In this case, each of the m n-by-p matrices in A is multiplied by B.
  Another example would be

  A[m-by-n-by-p] * D[m-p-by-q] = E[m-by-n-by-q]

  Now, each of the m n-by-p matrices in A is multiplied by the corresponding
  p-by-q matrix in B

2) be able to do the sort of thing Yorick does where, for example,
  A[n-by-m-by-p] * B[p-by-q-by-r] = C[n-by-m-by-q-by-r]

  This looks more like (2-dimensional) matrix mutiplication.

A simple way of achieving both would be to define a .dot(B,n) method
such that A.dot(B,n) would compute the dot product of the rows of A by
the corresponding rows of B (or any other dimension n, for that
matter), reducing the number of dimensions by one like add.reduce does
for a single array. The arrays should be conformable as for addition
and multiplication.

Essentially, A.dot(B,n) = add.reduce(A*B,n), but memory for the
intermediate result A*B would not need to be allocated (a saving when
casting of dimensions occurs, as in the example below).

For two-dimensional arrays, A.matrixMultiply(B) could be written
something like A[...,NewAxis].dot(B,-2).  The automatic casting of the
dimension introduced into A to match that of B would produce the
desired results. Note that add.reduce(A[...,NewAxis]*B) would imply
storing a 3-dimensional array. This would generalize well for any
number of dimensions, I believe. 

For example, case 1) above could be written

C = A[...,NewAxis].dot(B,-2)
E = A[...,NewAxis].dot(D[...,NewAxis,:,:],-2)

and case 2) could be written

A[...,NewAxis,NewAxis].dot(B,-3)

How do others feel about this issue?

Carlos Fonseca

P.S.: Why doesn't A[RubberIndex,NewAxis,NewAxis] work? 

	>>> a
	 0  1  2  3  4
	 5  6  7  8  9
	10 11 12 13 14
	>>> a[RubberIndex,NewAxis]
	Traceback (innermost last):
	  File "<stdin>", line 1, in ?
	ArrayError: each subindex must be either a slice, an integer,
	RubberIndex, or None
        ^^^^^^^^^^^ !??

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

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