np.dot 3x3 with N 1x3 arrays

Question

I have an ndarray of N 1x3 arrays I'd like to perform dot multiplication with a 3x3 matrix. I can't seem to figure out an efficient way to do this, as all the multi_dot and tensordot, etc methods seem to recursively sum or multiply the results of each operation. I simply want to apply a dot multiply the same way you can apply a scalar. I can do this with a for loop or list comprehension but it is much too slow for my application.

N = np.asarray([[1, 2, 3], [4, 5, 6], [7, 8, 9], ...])
m = np.asarray([[10, 20, 30], [40, 50, 60], [70, 80, 90]])

I'd like to perform something such as this but without any python loops:

np.asarray([np.dot(m, a) for a in N])

so that it simply returns [m * N[0], m * N[1], m * N[2], ...]

What's the most efficient way to do this? And is there a way to do this so that if N is just a single 1x3 matrix, it will just output the same as np.dot(m, N)?

Answer 1

Try This:

import numpy as np
N = np.asarray([[1, 2, 3], [4, 5, 6], [7, 8, 9], [1, 2, 3], [4, 5, 6]])
m = np.asarray([[10, 20, 30], [40, 50, 60], [70, 80, 90]])
re0 = np.asarray([np.dot(m, a) for a in N])  # original
re1 = np.dot(m, N.T).T  # efficient
print("result0:\n{}".format(re0))
print("result1:\n{}".format(re1))
print("Is result0 == result1? {}".format(np.array_equal(re0, re1)))

Output:

result0:
[[ 140  320  500]
 [ 320  770 1220]
 [ 500 1220 1940]
 [ 140  320  500]
 [ 320  770 1220]]
result1:
[[ 140  320  500]
 [ 320  770 1220]
 [ 500 1220 1940]
 [ 140  320  500]
 [ 320  770 1220]]
Is result0 == result1? True

Time cost:

import timeit
setup = '''
import numpy as np
N = np.random.random((1, 3))
m = np.asarray([[10, 20, 30], [40, 50, 60], [70, 80, 790]])
'''

>> timeit.timeit("np.asarray([np.dot(m, a) for a in N])", setup=setup, number=100000)
0.295798063278
>> timeit.timeit("np.dot(m, N.T).T", setup=setup, number=100000)
0.10135102272
# N = np.random.random((10, 3))
>> timeit.timeit("np.asarray([np.dot(m, a) for a in N])", setup=setup, number=100000)
1.7417007659969386
>> timeit.timeit("np.dot(m, N.T).T", setup=setup, number=100000)
0.1587108800013084
# N = np.random.random((100, 3))
>> timeit.timeit("np.asarray([np.dot(m, a) for a in N])", setup=setup, number=100000)
11.6454949379
>> timeit.timeit("np.dot(m, N.T).T", setup=setup, number=100000)
0.180465936661

Answer 2

First, regarding your last question. There's a difference between a (3,) N and (1,3):

In [171]: np.dot(m,[1,2,3])
Out[171]: array([140, 320, 500])       # (3,) result

In [172]: np.dot(m,[[1,2,3]])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-172-e8006b318a32> in <module>()
----> 1 np.dot(m,[[1,2,3]])

ValueError: shapes (3,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)

Your iterative version produces a (1,3) result:

In [174]: np.array([np.dot(m,a) for a in [[1,2,3]]])
Out[174]: array([[140, 320, 500]])

Make N a (4,3) array (this helps keep the first dim of N distinct):

In [176]: N = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10,11,12]])
In [177]: N.shape
Out[177]: (4, 3)
In [178]: np.array([np.dot(m,a) for a in N])
Out[178]: 
array([[ 140,  320,  500],
       [ 320,  770, 1220],
       [ 500, 1220, 1940],
       [ 680, 1670, 2660]])

Result is (4,3).

A simple dot doesn't work (same as in the (1,3) case):

In [179]: np.dot(m,N)
...
ValueError: shapes (3,3) and (4,3) not aligned: 3 (dim 1) != 4 (dim 0)
In [180]: np.dot(m,N.T)     # (3,3) dot with (3,4) -> (3,4)
Out[180]: 
array([[ 140,  320,  500,  680],
       [ 320,  770, 1220, 1670],
       [ 500, 1220, 1940, 2660]])

So this needs another transpose to match your iterative result.

The explicit indices of einsum can also take care of these transpose:

In [181]: np.einsum('ij,kj->ki',m,N)
Out[181]: 
array([[ 140,  320,  500],
       [ 320,  770, 1220],
       [ 500, 1220, 1940],
       [ 680, 1670, 2660]])

Also works with the (1,3) case (but not with the (3,) case):

In [182]: np.einsum('ij,kj->ki',m,[[1,2,3]])
Out[182]: array([[140, 320, 500]])

matmul, @ is also designed to calculate repeated dots - if the inputs are 3d (or broadcastable to that):

In [184]: (m@N[:,:,None]).shape
Out[184]: (4, 3, 1)
In [185]: (m@N[:,:,None])[:,:,0]     # to squeeze out that last dimension
Out[185]: 
array([[ 140,  320,  500],
       [ 320,  770, 1220],
       [ 500, 1220, 1940],
       [ 680, 1670, 2660]])

dot and matmul describe what happens with 1, 2 and 3d inputs. It can take some time, and experimentation, to get a feel for what is happening. The basic rule is last of A with 2nd to the last of B.

Your N is actually (n,3), n (3,) arrays. Here's what 4 (1,3) arrays looks like:

In [186]: N1 = N[:,None,:]
In [187]: N1.shape
Out[187]: (4, 1, 3)
In [188]: N1
Out[188]: 
array([[[ 1,  2,  3]],

       [[ 4,  5,  6]],

       [[ 7,  8,  9]],

       [[10, 11, 12]]])

and the dot as before (4,1,3) dot (3,3).T -> (4,1,3) -> (4,3)

In [190]: N1.dot(m.T).squeeze()
Out[190]: 
array([[ 140,  320,  500],
       [ 320,  770, 1220],
       [ 500, 1220, 1940],
       [ 680, 1670, 2660]])

and n of those:

In [191]: np.array([np.dot(a,m.T).squeeze() for a in N1])
Out[191]: 
array([[ 140,  320,  500],
       [ 320,  770, 1220],
       [ 500, 1220, 1940],
       [ 680, 1670, 2660]])