In Einstein summation notation, what numpy.dot() does now is:
c_riqk = a_rij * b_qjk

And you want:
c_[r]ik = a_[r]ij * b_[r]jk

where the brackets indicate a 'no summation' index.
Writing the ESN makes it clearer to me anyway.  :-)

--bb



On 5/18/06, Pau Gargallo <pau.gargallo@gmail.com> wrote:
> I'm afraid I really don't understand the operation that you want.

I think that the operation Angus wants is the following (at least I
would like that one ;-)

if you have two 2darrays of shapes:
     a.shape = (n,k)
    b.shape = (k,m)
you get:
    dot( a,b ).shape == (n,m)

Now, if you have higher dimensional arrays  (kind of "arrays of matrices")
    a.shape = I+(n,k)
    b.shape = J+(k,m)
where I and J are tuples, you get
    dot( a,b ).shape == I+J+(n,m)
    dot( a,b )[ i,j ] == dot( a[i],b[j] )          #i,j represent tuples
That is the current behaviour, it computes the matrix product between
every possible pair.
For me that is similar to 'outer' but with matrix product.

But sometimes it would be useful (at least for me) to have:
     a.shape = I+(n,k)
     b.shape = I+(k,m)
and to get only:
     dot2( a,b ).shape == I+(n,m)
     dot2( a,b )[i] == dot2( a[i], b[i] )

This would be a natural extension of the scalar product (a*b)[i] == a[i]*b[i]
If dot2 was a kind of ufunc, this will be the expected behaviour,
while the current dot's behaviour will be obtained by dot2.outer(a,b).

Does this make any sense?

Cheers,
pau