numpy tensordot related issue

2
votes

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?

3
pythonnumpymatrix-multiplication

3 Answers

3
votes

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]])
1
votes

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.

0
votes

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.