0

I have two numpy arrays A and B, both with the dimension [2,2,n], where n is a very large number. I want to matrix multiply A and B in the first two dimensions to get C, i.e. C=AB, where C has the dimension [2,2,n].

The simplest way to accomplish this is by using for loop, i.e.

for i in range(n):
    C[:,:,i] = np.matmul(A[:,:,i],B[:,:,i])

However, this is inefficient since n is very large. What's the most efficient way to do this with numpy?

  • https://stackoverflow.com/questions/40034993/how-to-get-element-wise-matrix-multiplication-hadamard-product-in-numpy – Sreeram TP Feb 23 '21 at 08:52
  • @SreeramTP That is only for 2D arrays right ? – joostblack Feb 23 '21 at 09:03
  • If that `i/n` dimension was first, you could just do `C =A@B`. In `numpy`, the leading dimension is outer most, so it's most natural to treat it as the "batch" dimension in `matmul`. – hpaulj Feb 23 '21 at 16:13

2 Answers2

0

You can do the following:

new_array = np.einsum('ijk,jlk->ilk', A, B)
joostblack
  • 2,465
  • 5
  • 14
0

What you want is the the default array multiplication in Numpy

In [22]: a = np.arange(8).reshape((2,2,2))+1 ; a[:,:,0], a[:,:,1]
Out[22]: 
(array([[1, 3],
        [5, 7]]),
 array([[2, 4],
        [6, 8]]))

In [23]: aa = a*a ; aa[:,:,0], aa[:,:,1]
Out[23]: 
(array([[ 1,  9],
        [25, 49]]),
 array([[ 4, 16],
        [36, 64]]))

Notice that I emphasized array because Numpy's arrays look like matrices but are indeed Numpy's ndarrays.


Post Scriptum

I guess that what you really want are matricesarrays with shape (n,2,2), so that you can address individual 2×2 matrices using a single index, e.g.,

In [27]: n = 3
    ...: a = np.arange(n*2*2)+1 ; a_22n, a_n22 = a.reshape((2,2,n)), a.reshape((n,2,2))
    ...: print(a_22n[0])
    ...: print(a_n22[0])
[[1 2 3]
 [4 5 6]]
[[1 2]
 [3 4]]

Post Post Scriptum

Re semantically correct:

In [13]: import numpy as np
    ...: n = 3
    ...: a = np.arange(2*2*n).reshape((2,2,n))+1
    ...: p = lambda t,a,n:print(t,*(a[:,:,i]for i in range(n)),sep=',\n')
    ...: p('Original array', a, n)
    ...: p('Using `einsum("ijk,jlk->ilk", ...)`', np.einsum('ijk,jlk->ilk', a, a), n)
    ...: p('Using standard multiplication', a*a, n)
Original array,
[[ 1  4]
 [ 7 10]],
[[ 2  5]
 [ 8 11]],
[[ 3  6]
 [ 9 12]]
Using `einsum("ijk,jlk->ilk", ...)`,
[[ 29  44]
 [ 77 128]],
[[ 44  65]
 [104 161]],
[[ 63  90]
 [135 198]]
Using standard multiplication,
[[  1  16]
 [ 49 100]],
[[  4  25]
 [ 64 121]],
[[  9  36]
 [ 81 144]]
gboffi
  • 22,939
  • 8
  • 54
  • 85