numpy tensordot related issue

Question

I have a specific issue with multiplying matrices in numpy. Here is an example:

P=np.arange(30).reshape((-1,3))
array([[ 0,  1,  2],
   [ 3,  4,  5],
   [ 6,  7,  8],
   [ 9, 10, 11],
   [12, 13, 14],
   [15, 16, 17],
   [18, 19, 20],
   [21, 22, 23],
   [24, 25, 26],
   [27, 28, 29]])

I want to multiply each row by its transpose in order to obtain a 3x3 matrix for each row, for example for the first row:

P[0]*P[0][:,np.newaxis]
array([[0, 0, 0],
   [0, 1, 2],
   [0, 2, 4]])

and store the result in a 3-d matrix M:

M=np.zeros((10,3,3))
for i in range(10):
    M[i] = P[i]*P[i][:,np.newaxis]

I think there might be a way to do this without looping, maybe with tensor-dot, but cannot find it.

Does someone have an idea?

score 3 · Accepted Answer · answered Mar 08 '11 at 08:57

3

It's just simple like this:

In []: P= arange(30).reshape(-1, 3)
In []: P[:, :, None]* P[:, None, :]
Out[]:
array([[[  0,   0,   0],
        [  0,   1,   2],
        [  0,   2,   4]],
       [[  9,  12,  15],
        [ 12,  16,  20],
        [ 15,  20,  25]],
       [[ 36,  42,  48],
        [ 42,  49,  56],
        [ 48,  56,  64]],
       #...   
       [[729, 756, 783],
        [756, 784, 812],
        [783, 812, 841]]])    
In []: P[1]* P[1][:, None]
Out[]:
array([[ 9, 12, 15],
       [12, 16, 20],
       [15, 20, 25]])

answered Mar 08 '11 at 08:57

eat

7,440
1
19
27

1

great, exactly what I was looking for, is there an easy methodology to understand which index I should expand with None? – Andrea Zonca Mar 08 '11 at 16:58
@Andrea Z: For such methodology; I'll first try to figure out how the `numpy` broadcasting really works. Then just realizing you are looking for a `shape= (10, 3, 3)´, it would be more or less 'natural' to figure out that `(P[:, :, None]* P[:, None, :]).shape== (10, 3, 3)`. Perhaps not a best description, but the main point I stand for is; familiarize your self how broadcasting works. Now feel free to ask an other question on 'the innards of broadcasting'! Thanks – eat Mar 08 '11 at 18:24

score 1 · Answer 2 · edited May 23 '17 at 11:55

1

Since I like stride_tricks, that's what I'd use. I'm sure there are other ways.

Change the stride and shape of the array so that you expand it into 3D. You could easily do the same thing with the "transposed" version of P, but here I just reshape it and let the broadcasting rules stretch it into the other dimension.

P=np.arange(30).reshape((-1,3))
astd = numpy.lib.stride_tricks.as_strided
its = P.itemsize
M = astd(P,(10,3,3),(its*3,its,0))*P.reshape((10,1,3))

I'm going to add a reference to this post because it is a nice detailed explanation of stride_tricks.as_strided.

edited May 23 '17 at 11:55

Community

1
1

answered Mar 08 '11 at 05:50

Paul

42,322
15
106
123

+1 stride_tricks is really wonderful, but so complex to grasp! – Andrea Zonca Mar 08 '11 at 17:00

score 0 · Answer 3 · answered Mar 08 '11 at 09:50

0

This partially solves the problem using tensordot(),

from numpy import arange,tensordot

P = arange(30).reshape((-1,3))

i = 3

T = tensordot(P,P,0)[:,:,i,:]

print T[i]
print tensordot(P[i],P[i],0)

T contains all the products you want (and more), it's just a question of extracting them.

answered Mar 08 '11 at 09:50

lafras

8,712
4
29
28

Although OP speculated the possibility to utilize `tensordot(.)`, nothing in his code actually requires it. Your solution seems to be very expensive both execution time and memory consumption wise. Have you made any benchmarks to compare yours solution with the other ones? Thanks – eat Mar 08 '11 at 14:45
@eat: True, it is expensive memory wise. In this example, with a 10 by 3 array, execution time is approximately double using `tensordot()`. It gets progressively worse as the array size increases. – lafras Mar 08 '11 at 20:21
1

I didn't make the point explicitly, but tensordot is probably _not_ the way to do it, for a completely different reason. If you view `P` as a rank-2 tensor then only three options exist for the product of `P` with itself, 1) either all the indexes cancel leaving you with a rank-0 tensor (a scalar), 2) 1 set of indexes cancels and you are left with a rank-2 tensor (a matrix), or 3) none of them cancel and you're left with a rank-4 tensor. The OP wanted a rank-3 tensor as an answer, which `tensordot()` won't give you, unless you slice the rank-4 product. – lafras Mar 08 '11 at 20:21

numpy tensordot related issue

3 Answers3