Python - Matrix outer product

Question

Given two matrices

A: m * r
B: n * r

I want to generate another matrix C: m * n, with each entry C_ij being a matrix calculated by the outer product of A_i and B_j.

For example,

A: [[1, 2],
    [3, 4]]

B: [[3, 1],
    [1, 2]]

gives

C: [[[3, 1],  [[1 ,2],
     [6, 2]],  [2 ,4]],
     [9, 3],  [[3, 6],
     [12,4]],  [4, 8]]]

I can do it using for loops, like

    for i in range (A.shape(0)):
      for j in range (B.shape(0)):
         C_ij = np.outer(A_i, B_j)

I wonder If there is a vectorised way of doing this calculation to speed it up?

Bomber

Do you want a 4D, (m, n, r, r)-shape array, or do you want a 2D, (m, n)-shape array of object dtype where each element is another array? I would strongly recommend the first option, but your description sounds closer to the second. – Semiannual

Sorry for the confusion, but I prefer the first one, a 4D (m, n, r, r)-shape array. – Bomber

Answer 1

temp = numpy.multiply.outer(A, B)
C = numpy.swapaxes(temp, 1, 2)

NumPy ufuncs, such as multiply, have an outer method that almost does what you want. The following:

temp = numpy.multiply.outer(A, B)

produces a result such that temp[a, b, c, d] == A[a, b] * B[c, d]. You want C[a, b, c, d] == A[a, c] * B[b, d]. The swapaxes call rearranges temp to put it in the order you want.

Semiannual

Answer 2

The Einstein notation expresses this problem nicely

In [85]: np.einsum('ac,bd->abcd',A,B)
Out[85]: 
array([[[[ 3,  1],
         [ 6,  2]],

        [[ 1,  2],
         [ 2,  4]]],


       [[[ 9,  3],
         [12,  4]],

        [[ 3,  6],
         [ 4,  8]]]])

Nightrider

Man, I ought to learn that notation. Every time someone posts an answer using it, it's always way shorter than what I come up with. Probably more understandable too, if you know the Einstein summation convention. – Semiannual

Answer 3

Simple Solution with Numpy Array Broadcasting

Since, you want C_ij = A_i * B_j, this can be achieved simply by numpy broadcasting on element-wise-product of column-vector-A and row-vector-B, as shown below:

# import numpy as np
# A = [[1, 2], [3, 4]]
# B = [[3, 1], [1, 2]]
A, B = np.array(A), np.array(B)
C = A.reshape(-1,1) * B.reshape(1,-1)
# same as: 
# C = np.einsum('i,j->ij', A.flatten(), B.flatten())
print(C)

Output:

array([[ 3,  1,  1,  2],
       [ 6,  2,  2,  4],
       [ 9,  3,  3,  6],
       [12,  4,  4,  8]])

You could then get your desired four sub-matrices by using numpy.dsplit() or numpy.array_split() as follows:

np.dsplit(C.reshape(2, 2, 4), 2)
# same as:
# np.array_split(C.reshape(2,2,4), 2, axis=2)

Output:

[array([[[ 3,  1],
         [ 6,  2]],

        [[ 9,  3],
         [12,  4]]]), 
array([[[1, 2],
         [2, 4]],

        [[3, 6],
         [4, 8]]])]

Blackhead

Answer 4

You can use

C = numpy.tensordot(A, B, axes=0)

tensordot does exactly what you want. The axes parameter is used to in addition perform sums over certain axes (for >2d tensors, the default value is 2 and it will sum away 2 axes of each of the arrays), but by setting it to 0 it simply doesn't reduce the dimensions and keeps the whole outer product.

Emulsion

Answer 5

Use numpy;

In [1]: import numpy as np

In [2]: A = np.array([[1, 2], [3, 4]])

In [3]: B = np.array([[3, 1], [1, 2]])

In [4]: C = np.outer(A, B)

In [5]: C
Out[5]: 
array([[ 3,  1,  1,  2],
       [ 6,  2,  2,  4],
       [ 9,  3,  3,  6],
       [12,  4,  4,  8]])

Once you have the desired result, you can use numpy.reshape() to mold it in almost any shape you want;

In [6]: C.reshape([4,2,2])
Out[6]: 
array([[[ 3,  1],
        [ 1,  2]],

       [[ 6,  2],
        [ 2,  4]],

       [[ 9,  3],
        [ 3,  6]],

       [[12,  4],
        [ 4,  8]]])

Lumbard

I'm kind of surprised this is so visually close to what the OP wants, but despite the resemblance, it's still not quite right. If you visualize the desired result as a big 2D grid with little 2D grids in the the cells, this is what you'd get by gluing all the little grids together at the edges to make one grid. – Semiannual

Yes, agree with you. Perhaps first do np.outer(A, B) and then divide it into smaller grids? – Bomber

But this will give a different result. I wanted the inner matrices to be [[3, 1], [6, 2]], [[1, 2], [2, 4]], [[9, 3], [12, 4]] and [[3, 6], [4, 8]]. – Bomber