I am having some problems relating to the current function dot
(identical to matrixmultiply, though I haven't seen the equivalence in
any documentation). Here is what the docs say:
dot(a,b)
Will return the dot product of a and b. This is equivalent to matrix
multiply for 2d arrays (without the transpose). Somebody who does
more linear algebra really needs to do this function right some day!
Or the builtin doc string:
>>> print Numeric.dot.__doc__
dot(a,b) returns matrix-multiplication between a and b. The product-sum
is over the last dimension of a and the second-to-last dimension of b.
First, this is misleading. It seems to me to indicate that b must
have rank at least 2, which experiments indicate is not necessary.
Instead, the rule appears to be to use the only axis of b if b has
rank 1, and otherwise to use the second-to-last one.
Frankly, I think this convention is ill motivated, hard to remember,
and even harder to justify. As a mathematician, I can see only one
reasonable default choice: One should sum over the last index of a,
and the first index of b. Using the Einstein summation convention
[*], that would mean that
dot(a,b)[j,k,...,m,n,...] = a[j,k,...,i] * b[i,m,n,...]
[*] that is, summing over repeated indices -- i in this example
This would of course yield the current behaviour in the important
cases where the rank of b is 1 or 2.
But we could do better than this: Why not leave the choice up to the
user? We could allow an optional third parameter which should be a
pair of indices, indicating the axes to be summed over. The default
value of this parameter would be (-1, 0). Returning to my example
above, the user could then easily compute, for example,
dot(a,b,(1,2))[j,k,...,m,n,...] = a[j,i,k,...] * b[m,n,i,...]
while the current behaviour of dot would correspond to the new
behaviour of dot(a,b,(-1,-2)) whenever b has rank at least 2.
Actually, there is probably a lot of code out there that uses the
current behaviour of dot. So I would propose leaving dot untouched,
and introducing inner instead, with the behaviour I outlined above.
We could even allow any number of pairs of axes to be summed over, for
example
inner(a,b,(1,2),(2,0))[k,l,...,m,n,...] = a[k,i,j,l,..] * b[j,m,i,n,...]
With this notation, one can for example write the Hilbert-Schmidt
inner product of two real 2x2 matrices (the sum of a[i,j]b[j,i] over
all i and j) as inner(a,b,(0,1),(1,0)).
If my proposal is accepted, the documentation should probably declare
dot (and its alias matrixmultiply?) as deprecated and due to disappear
in the future, with a pointer to its replacement inner. In the
meantime, dot could in fact be replaced by a simple wrapper to inner:
def dot(a,b):
if len(b.shape) > 1:
return inner(a,b,(-1,-2)
else:
return inner(a,b)
(with the proper embellishments to allow it to be used with python
sequences, of course).
- Harald
