outerproduct

a1 = np.array([[1, 2, 3], [1, 2, 3]])
a2 = np.array([[1, 2], [1, 2]])
a3 = np.array([[1, 2, 3, 4], [1, 2, 3, 4]])
 
np.einsum('bi,bj,bk->bijk', a1, a2, a3).reshape(2, -1).shape