I am sorry, but I don't agree with you.
The matrix class seems to be there to help people with numerical linear
algebra algorithms "a la matlab". Is this an incorrect assumption?
The behavior I am asking can be restricted to the matrix class objects
only. However as dot already makes a scalar out the multiplication of
two rank-1 arrays (in which case it computes the inner product), I
thought that this behavior could be extended to matrix objects. Another
possibility (without touching dot) would be to change the functions
__mul__, __rmul__, and __imul__ in the matrix class to deal with 1x1
matrices as scalars. That would be OK to me too.
The current behavior is odd when writing numerical algorithms. Take a
look at these two examples:
1) Computing the norm:
The natural way:
def norm(x): return sqrt(x.T*x)
But in this case you can't write norm(v)*v!
Then, you have then to define norm as
def norm(x): return sqrt(x.T*x).item()
2) The conjugate gradient iteration:
alpha = (r.T*r) / (p.T*A*p)
x += alpha*p
rold = rnew
r += alpha*(A*p)
beta = (r.T*r) / (rold.T*rold)
p = -r + beta*p
This returns an error on second, forth, and the last equation! In all
these cases, you'll need a call to the "item" method.
Finally, the example you gave
(x*x.T)*(x.T*x)
does not make sense as matrix multiplication formula unless you
interpret the first part as a scalar. You seem to be happy that the
compiler got "the mistake". But in any Mathematics book the formula
above would be considered OK as the first part would be interpreted as
an inner product (and the second an outer product). In my opinion, it is
not a mistake in a linear algebra setting.
My 2 cents,
Paulo
--
Paulo José da Silva e Silva
Professor Assistente do Dep. de Ciência da Computação
(Assistant Professor of the Computer Science Dept.)
Universidade de São Paulo - Brazil
e-mail: rsilva at ime.usp.br Web: http://www.ime.usp.br/~rsilva
Teoria é o que não entendemos o (Theory is something we don't)
suficiente para chamar de prática. (understand well enough to call
practice)