Add Multiple Of A Matrix Without Build A New One

Say I have two matrices B and M and I want to execute the following statement: B += 3*M I execute this instruction repeatedly so I don't want to build each time the matrix 3*M (3

Solution 1:

You can use np.tensordot -

np.tensordot(As,Ms,axes=(0,0))

Or np.einsum -

np.einsum('i,ijk->jk',As,Ms)

Sample run -

In [41]: As = [2,5,6]

In [42]: Ms = [np.random.rand(2,3),np.random.rand(2,3),np.random.rand(2,3)]

In [43]: sum(a*M for a, M in zip(As, Ms))
Out[43]: 
array([[  6.79630284,   5.04212877,  10.76217631],
       [  4.91927651,   1.98115548,   6.13705742]])

In [44]: np.tensordot(As,Ms,axes=(0,0))
Out[44]: 
array([[  6.79630284,   5.04212877,  10.76217631],
       [  4.91927651,   1.98115548,   6.13705742]])

In [45]: np.einsum('i,ijk->jk',As,Ms)
Out[45]: 
array([[  6.79630284,   5.04212877,  10.76217631],
       [  4.91927651,   1.98115548,   6.13705742]])

Solution 2:

Another way you could do this, particularly if you favour readability, is to make use of broadcasting.

So you could make a 3D array from the 1D and 2D arrays and then sum over the appropriate axis:

>>> Ms = np.random.randn(4, 2, 3)   # 4 arrays of size 2x3
>>> As = np.random.randn(4)
>>> np.sum(As[:, np.newaxis, np.newaxis] * Ms)
array([[-1.40199248, -0.40337845, -0.69986566],
       [ 3.52724279,  0.19547118,  2.1485559 ]])
>>> sum(a*M for a, M in zip(As, Ms))
array([[-1.40199248, -0.40337845, -0.69986566],
       [ 3.52724279,  0.19547118,  2.1485559 ]])

However, it's worth noting that np.einsum and np.tensordot are usually much more efficient:

>>> %timeit np.sum(As[:, np.newaxis, np.newaxis] * Ms, axis=0)
The slowest run took 7.38 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 8.58 µs per loop
>>> %timeit np.einsum('i,ijk->jk', As, Ms)
The slowest run took 19.16 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 2.44 µs per loop

And this is also true for larger numbers:

>>>Ms = np.random.randn(100, 200, 300)>>>As = np.random.randn(100)>>>%timeit np.einsum('i,ijk->jk', As, Ms)
100 loops, best of 3: 5.03 ms per loop
>>>%timeit np.sum(As[:, np.newaxis, np.newaxis] * Ms, axis=0)
100 loops, best of 3: 14.8 ms per loop
>>>%timeit np.tensordot(As,Ms,axes=(0,0))
100 loops, best of 3: 2.79 ms per loop

So np.tensordot works best in this case.

The only good reason to use np.sum and broadcasting is to make the code a little more readable (helps when you have small matrices).