Operators for matrix: current choices (Was Matlab vs Python ...)

[Tim Hochberg]

hzhu at localhost.localdomain (Huaiyu Zhu) writes:

[SNIP some stuff about how proposed notation is not consisten with
existing NumPy behaviour]

> I would have prefered compatibility with NumPy, except for two reasons:
> 1. MatPy was started because there is a real need for using infix operator
>    for matrix multiplication and the only way available was to define a new
>    class.  So for MatPy users, at least, * already means matrixwise.
> 2. The current list add is not componentwise.

This is fair enough. You could add that the existing list replication
(*) is not componentwise either. However, with this incompatability in
mind I think "just tell the newbies that . means componentwise" is a
bit of an overstatement. (I'm afraid I don't have the original in
front of me, so I've probably butchered that quote, sorry).

> I don't know of a way that satisfies every compatibility demand.  Because
> the supply of binary operators has been very limited, each package assigns
> them to their most deserving ones instead of according to a better overall
> plan.  Had there been more this could have been better coordinated.
> 
> >Personally I think that this an overlarge crop of operators to be
> >adding. Since I prefer something like:
> >
> >A.inv * B 
> >
> >to 
> >
> >A \ B
> 
> They are equal in math, but not in practice.  The former takes O(n^3)
> computation while the latter is O(n^2).  The former is also less accurate
> for near singular A because of rounding error.

Actually it's easy enough to make them equal in practice as well as in
math. A.inv doesn't have to be (in fact shouldn't be) the actual
inverse matrix, rather it starts it's life as a class containing a
reference to A. As a simple example:

class InverseMatrix:

   def __init__(self, matrix):
      self._matrix = matrix

   def __mul__(self, other):
      return LinearAlgebra.solve_linear_equations(self._matrix, other)

where A.inv is initialize to InverseMatrix(A).

Then A.inv * B returns LinearAlgebra.solve_linear_equation(A, B) and
is order N**2 just as A\B.

Things would have to more complicated in real life of course. If
someone tried to look at or modify A.inv, you'd probably want to cache
the modified result. Some thought would be required to get this right,
but I believe it's doable.

> The other operators are also useful. For example, the + and .+ do different
> things.  If you add a row and column vector with + it will complain about
> shape mismatch, as defined in linear algebra. With .+ it will give a matrix
> properly extended to both rows and columns.  This has been asked for in this
> thread. I've seen it been called prolongation, continuation or cross
> extension, among others.  In fact, in MatPy, the .+ .- .* ./ do exactly what
> + - * / do in NumPy.

I'm afraid I don't find this usage all that compelling, but that, of
course, is just my opinion.

-tim