Re: [Numpy-discussion] big picture? One proposal

[Pearu P. <pe...@ce...>]

Hi,

On Fri, 8 Mar 2002, Konrad Hinsen wrote:

<removed relevant comments>

> > Numeric community is united on this, I think Guido would be receptive.
> > We might suggest a particular operator symbol or pair (triple) but 
> 
> Actually I feel quite safe: there might be a majority for another
> operator, but I don't expect we'd ever agree on a symbol :-)

Thanks Konrad for this excellent point. Seeing all these proposals about
solving our matrix-array issues makes also me feel safer with the current
situation. My general point is that a _good_ solution is simple,
if a solution is not simple, then it is probably a bad solution.

I find separating array and matrix instances (in a sense of raising
exception when doing <matrix> <op> <array>) not a very simple solution:
  New concepts are introduced that actually do not solve the simplicity
problem of representing matrix operations. As I see it, they
only introduce restrictions and the main assumption behind the rationale
is that "users are dumb and they don't know what is best for them".
This is how I interpret the raised exception as behind the scenes matrix
and array are the same (in the sense of data representation). 

Let me remind at least to myself that the discussion started from a
proposal that aimed to write matrix multiplication of arrays

  Numeric.dot(a,b)

in somewhat simpler form. Current solution is to use Matrix class that
being just a wrapper of arrays, redefines __mul__ method; and after
one has defined
  a = Matrix.Matrix(a)
the matrix multiplication of arrays looks simple:
  a * b

Travis's proposal was to reduce the first step and have it inside an
expression in a short form:

  a.M * b

(there have been two implementation approaches proposed for this: (i) a.M
returns a Matrix instance, (ii) a.M returns the same array with a
temporarily set bit saying that the following operation is somehow
special).
To me, this looks like a safe solution. Though it is a hack, at least it
is simple and understandable in anywhere where it is used (having a * b
where b can be either matrix or array, it is not predictable from just
looking the code what the result will be -- not very pythonic indeed).

The main objection to this proposal seems to be that it deviates from a
good pythonic style (ie don't mess with attributes in this way).
I'd say that if python does not provide a good solution to our problem,
then we are entitled to deviate from a general style. After all, in doing
numerics the efficiency issue has a rather high weight. And a generally
good style of Python cannot always support that.


I guess I am missing here a goal to get Numeric or numarray joined
to Python. With this perspective the only efficient solution seems
to be introducing a new operator (or new operators). Few candidates have
been proposed:

 a ~* b       - BTW, to me this looks like dot(conjugate(a),b).
 a (*) b      - note the in situ version of it: a (*)= b
 a (**) b     - looks ugly enough? ;-)
Actually, why not

 a [*] b, a {*} b

for direct products of matrices (BTW (*) seems more appropriate here). So,
my ideal preference would be:

  a .* b   - element-wise multiplication of arrays, 2nd pref.: a * b
   a * b   - matrix multiplication of arrays, 2nd preference: a [*] b
  a (*) b  - direct matrix multiplication (also know as tensor product)
             of arrays
    a~     - conjugate of arrays
    a`     - transpose of arrays

This looks great but requires many new features to Python (new operators, 
the concept of right-hand unary operator).

I don't think that Python should introduce these new operators just
because of Numeric community. It is fine if they get used in other fields
as well that suffer from the lack of operators.


About unary operations: transpose and conjugate. BTW, in complex linear
algebra their composition is equally frequent operator. Let me propose the
following solution: to have

  a ** T  for Numeric.transpose(a)
  a ** H  for Numeric.transpose(Numeric.conjugate(a))

define

T = TransposeOp()
H = TransposeOp(conjugate=1)

where

class TransposeOp:
    def __init__(self, conjugate=0):
         self.conjugate = conjugate
    def __rpow__(self,arr):
         if self.conjugate:
             return Numeric.transpose(Numeric.conjugate(a))
         return Numeric.transpose(arr)

Looks Pythonic to me;-)

Regards,
	Pearu