[Numpy-discussion] Re: Dot product threading?

[Robert Kern]

Angus McMorland wrote:
> Robert Kern wrote:
> 
>>Angus McMorland wrote:
>>
>>>Is there a way to specify which dimensions I want dot to work over?
>>
>>Use swapaxes() on the arrays to put the desired axes in the right places.
> 
> Thanks for your reply, Robert. I've explored a bit further, and have
> made sense of what's going on, to some extent, but have further questions.
> 
> My interpretation of the dot docstring, is that the shapes I need are:
>     a.shape == (2,3,1) and b.shape == (2,1,3)
> so that the sum is over the 1s, giving result.shape == (2,3,3)

I'm not sure why you think you should get that resulting shape. Yes, it will
"sum" over the 1s (in this case there is only one element in those axes so there
is nothing really to sum). What exactly are the semantics of the operation that
you want? I can't tell from just the input and output shapes.

> but:
>     In [85]:ma = array([[[4],[5],[6]],[[7],[8],[9]]]) #shape = (2, 3, 1)
>     In [86]:mb = array([[[4,5,6]],[[7,8,9]]])         #shape = (2, 1, 3)
> so
>     In [87]:res = dot(ma,mb).shape
>     In [88]:res.shape
>     Out[88]:(2, 3, 2, 3)
> such that
>     res[i,:,j,:] == dot(ma[i,:,:], mb[j,:,:])
> which means that I can take the results I want out of res by slicing
> (somehow) res[0,:,0,:] and res[1,:,1,:] out.
> 
> Is there an easier way, which would make dot only calculate the dot
> products for the cases where i==j (which is what I call threading over
> the first dimension)?

I'm afraid I really don't understand the operation that you want.

> Since the docstring makes no mention of what happens over other
> dimensions, should that be added, or is this the conventional numpy
> behaviour that I need to get used to?

It's fairly conventional for operations that reduce values along an axis to a
single value. The remaining axes are left untouched. E.g.

In [1]: from numpy import *

In [2]: a = random.randint(0, 10, size=(3,4,5))

In [3]: s1 = sum(a, axis=1)

In [4]: a.shape
Out[4]: (3, 4, 5)

In [5]: s1.shape
Out[5]: (3, 5)

In [6]: for i in range(3):
   ...:     for j in range(5):
   ...:         print i, j, (sum(a[i,:,j]) == s1[i,j]).all()
   ...:
   ...:
0 0 True
0 1 True
0 2 True
0 3 True
0 4 True
1 0 True
1 1 True
1 2 True
1 3 True
1 4 True
2 0 True
2 1 True
2 2 True
2 3 True
2 4 True

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco