Discussion: new operators for numerical computation

[Tim Hochberg]

Huaiyu Zhu <huaiyu_zhu at yahoo.com> writes:

> Now that Guido has given the go-ahead for adding new operators for numerical
> computation, let us consider carefully how to allocate this new resource,
> taken into account Tim's guidelines.
> 
> Here's a list of issues that come to mind:
> 
> 1. What symbols to use.
> 
>    It seems that we can more or less concentrate on using either . or @
>    together with existing math operators.  Issues:
> 
>    (1) Use . or @

 Not '@'. It's too dense and looks to much like a letter. '.' is
 better, but is has the potential to cause confusion in situations
 such as 3.+M. This equation is ambiguous if the ".x" operators are
 elementwise, but not if they're matrixwise (I'm assuming that the
 elementwise operators obey the same broadcasting rules that NumPy's
 operators do today). 

>    (2) Which is matrixwise or elementwise

I can't see a good answer here. NumPy has been around for a long time
and I believe it has a fairly large user base which is a strong
argument in favor of the ".x" operators being matrixwise. However,
other operators within python (*,+) operate "objectwise" on sequences
(concatenation and replication). I think it would be suprising if
these operators were elementwise for matrices and objectwise for
other, similar, objects (sequences). This tends to favor ".x" being
elementwise.

>    (3) Compatibility: NumPy uses * for elementwise, MatPy uses * for
>        matrixwise.

Do you have an idea of what the relative sizes of the userbases of
these two packages are? Honest question -- I have no idea, but I've
been operating under the assumption that userbase of MatPy was
relatively small.

>    (4) Similarity with Matlab (but not really compatibility). 

I think it's a mistake to look just at Matlab here. If were going to
look to matrix languages for inspriation we should look to
several. Unfortunately, I'm no help here.


> 3. Their behavior on other objects:
> 
>    Even though they are supposed to be used for numerical objects only, what
>    if user writes [10,20].+[1,2]?
> 
>    - Exception: no __dotadd__ or __rdotadd__ defined
> 
>    - Behave as eadd and give [11,22]
> 
>    - Behave as current + and give [10,20,1,2]
> 
>    The last one is dead, people writing [].+[] either made a typo of
>    [].member+[] or want eadd.
> 
>    The second seems most sensible and least surprising, but it does need to
>    include list in the implementation of eadd.
> 
>    My understanding is that Tim wants the first one.  So let's leave it
>    there lest it drags other bits down.

Note that if ".x" ends up being matrixwise, things will become
horribly inconsistent if it is later decided to extend these operators
to other sequences. This IMO is probably the strongest argument for
making ".x" elementwise.

> 
> 
> 4. The magic words:
> 
>    This really depends on issue 3 above.  
> 
>    - If it is intended that these thing should never acquire a meaning
>      outside numerical computation, they should be called __mmul__,
>      __emul__, etc, because that would describe the intension well.
> 
>    - If other classes, like UserList, are allowed to use them, it's better
>      to use __dotmul__, __atmul__, etc, so the association with the symbols
>      are strong.
> 
>    - Someone mentioned __altmul__ for @*, meaning alternative mul.

Not ready to dive in here till some other stuff is resolved.

> 
> 5. Other bits and pieces:
> 
>    - How to deal with left division?
> 
>     (1) Use sol(a,b)
>     (2) use /@ or /.
>     (3) use %

'%' is already a valid elementwise operator and to have "%" and ".%"
mean completely different things is not good, so (3) is out. 

I think (2) only makes sense if the ".x" operators are matrixwise. It
seems just too weird to have "/" be right division "/." be left
division and "./" be elementwise. It's much more sensible to have "/"
be elementwise, "/." be left division, and "./" be right division for
instance.

If possible I would like to avoid (1). 'sol' would almost certainly be
in the a seperate module and this is potentially painful.

>    - How to deal with power?  It is already two characters
> 
>     (1) Use ** and .** or @**
>     (2) Steal ^ with .^ or @^ (not much chance, just mentioning here)

It's gotta be (1). Unless a limited set of operations is considered
(see below), then _maybe_ ".^".

>    - Do we want a copy operator := which is different from =?  It allows
>      efficient code.  The following code does not change A
>          B := A
>          B += C
>      This definitely will be usable for other classes. Call it __copy__?

This seems to be a general request and not particularly numeric
specific. And I don't think it adds much in the way of
functionality. I would drop this. 

>    - Since augmented assignment almost certainly will get in the core, we
>      need to consider combining them.  A few possible opinions
> 
>      (1) Just use prefix . or @ to mean e or m, use suffix = to mean in
>          place assignment, and all are well defined.
> 
>      (2) a .**= 3 looks too long.  Don't combine them.
> 
>      (3) allow .+=, etc but ban .**=.

Here I'd go with (1), consistency over beauty.


After rereading my answers to the above, I'm going to have to come
down on the side of ".x" being matrixwise, despite the clash with
sequence concatenation and replication. It removes the ambiguity of
the grammer, it makes more sense with "./" and "/." and it's consistent
with NumPy. 

I'm not yet convinced as yet that all of these operators are needed;
the +, - and .+, .- operators are very similar in behaviour. The main
reason they are present seems to be orthogonality. I'm sure this has
been brought up before, but does MatLab have both broadcasting and non
broadcasting versions of its operators? Is it worth considering adding
only the matrix operators .*, ./, /., and(?) .** and letting NumPy/
MatPy choose whether they wish to broadcast or not? It seems that the
difference between .* and * is different in character than the
proposed difference between .+ and +; in one it's between elementwise
and matrixwise operations, and in the second it's between broadcasting
and nonbroadcasting operations. Yes, you can think of the second in
terms of the first, but it still feels a bit forced to me.

Based on my ramblings I'll toss out a proposed operator set for the
pack to feed upon:

	Operator	Name
          .*             matmul
          ./             matdiv
          /.             matsol
          .** or(?) .^   matpow
	  *		 mul
	  /		 div
	  %,^,&,|,**     mod, xor, and, or, pow

With broadcasting / nonbroadcasting decided by the individual packages.

My 3 cents.

-tim