Numpy: broadcast translation on one common axis from two 2d arrays

Question

I am looking for a way to stepwise multiply two 2d arrays of form (a, b) and (b, c) respectively. The b axis, which has two arrays.

For example, an example of what I would like to broadcast (vectorize) is:

import numpy as np    

# some dummy data
A = np.empty((2, 3))
B = np.empty((3, 4))

# naive implementation
C = np.vstack(np.kron(A[:, i], B[i, :]) for i in [0, 1, 2])

# this should give (3, 2, 4)
C.shape

Does anyone know what to do here? Is there a better way?

Answer 1

credit for @hpaulj for definitions AandB
use np.outerandnp.stack

A = np.arange(6).reshape((2, 3))
B = np.arange(12).reshape((3, 4))

np.stack([np.outer(A[:, i], B[i, :]) for i in range(A.shape[1])])

[[[ 0  0  0  0]
  [ 0  3  6  9]]

 [[ 4  5  6  7]
  [16 20 24 28]]

 [[16 18 20 22]
  [40 45 50 55]]]

---

and get np.einsumin the right shape

np.einsum('ij, jk->jik', A, B)

[[[ 0  0  0  0]
  [ 0  3  6  9]]

 [[ 4  5  6  7]
  [16 20 24 28]]

 [[16 18 20 22]
  [40 45 50 55]]]

---

Broadcasting and transpose

(A[:, None] * B.T).transpose(2, 0, 1)

[[[ 0  0  0  0]
  [ 0  3  6  9]]

 [[ 4  5  6  7]
  [16 20 24 28]]

 [[16 18 20 22]
  [40 45 50 55]]]

---

the form (3, 2, 4)

time

Answer 2

In different test cases:

In [56]: A=np.arange(6).reshape((2,3))
In [57]: B=np.arange(12).reshape((3,4))
In [58]: np.vstack([np.kron(A[:,i],B[i,:]) for i in range(3)])
Out[58]: 
array([[ 0,  0,  0,  0,  0,  3,  6,  9],
       [ 4,  5,  6,  7, 16, 20, 24, 28],
       [16, 18, 20, 22, 40, 45, 50, 55]])

First try with `einsum, saving all 3 axes (no summation)

In [60]: np.einsum('ij,jk->ijk',A,B)
Out[60]: 
array([[[ 0,  0,  0,  0],
        [ 4,  5,  6,  7],
        [16, 18, 20, 22]],

       [[ 0,  3,  6,  9],
        [16, 20, 24, 28],
        [40, 45, 50, 55]]])

, .

, 2x4x3, 8,3 .

In [64]: np.einsum('ij,jk->ikj',A,B).reshape(8,3).T
Out[64]: 
array([[ 0,  0,  0,  0,  0,  3,  6,  9],
       [ 4,  5,  6,  7, 16, 20, 24, 28],
       [16, 18, 20, 22, 40, 45, 50, 55]])

,

In [68]: np.einsum('ij,jk->jik',A,B).reshape(3,8)
Out[68]: 
array([[ 0,  0,  0,  0,  0,  3,  6,  9],
       [ 4,  5,  6,  7, 16, 20, 24, 28],
       [16, 18, 20, 22, 40, 45, 50, 55]])

. A (2,3), B (3,4), (3,2,4) (3,8). = 2, j = 3, k = 4 = > jik.

, ,

a_ij * b_jk = c_jik

sum einsum, .