Multidimensional matrix multiplication in python

Question

I have matrix A of dimension 500x2000x30 and matrix B of dimension 30x5.

You can think that there are 500 instances of 2000x30 as matrix A is of dimension 500x2000x30.

I want to multiply each of 1x2000x30 from A with matrix B to obtain new matrix of size 1x2000x5.

i.e. A X B should give me a matrix of dimension 500x2000x5

Obviously looping 500 times through matrix A is a solution but is there an efficient way to achieve this?

Edit: Both A and B are numpy arrays

Answer 1

First of all any multidimensional array is a n times b split of a flatten array. The 2,3,2 dimension array, [ [ [A,B],[C,D],[E,F] ], [ [G,H],[I,J],[K,L] ] ] is simply divide A,B,C,D,E,F,G,H,I,J,K,L into 2 pieces , then 3 pieces, then 2 pieces again and there is your multidimensional array. Vectors are here two dimensional.

Secondly, two dimensional matrix multiplication is also fixed length vector multiplication with a combination. (n,m)*(m,l) dimensional matrix multiplication is actually term by term multiplication and sum of results of l different m vectors and n different m vectors. l times n combination of two different m vectors.

When you get that you simply can convert 20 dimension matrix to two dimension matrix with same vector size and multiply them and reshape it back.

Quick Example:

import numpy as np

a=np.random.random(120).reshape(2,5,3,4)
b=np.random.random(120).reshape(5,3,4,2)

br=np.rollaxis(b,3,2).reshape(30,4).T  
# converted all combinations to a series of 4 dimension vectors here 

test1=np.matmul(a.reshape(30,4),br).reshape(2,5,3,5,3,2)
test2=np.dot(a,b)

np.all(test1==test2)

returns

array(True)

lets do it basically for all combinations (dimensions)

sa=a.shape
sb=b.shape
res=np.ndarray([sa[0],sa[1],sb[0],sb[1],sa[2],sb[3]])
for i in range(a.shape[0]):
  for j in range(a.shape[1]):
    for k in range(b.shape[0]):
      for n in range(b.shape[1]):
        x=np.matmul(a[i,j],b[k,n])
        np.append(res,x)
print(np.any(res==np.dot(a,b).reshape(2,5,5,3,3,2)))

returns

True

USE CASE:

Suppose we have a case such as for 10 different materials and 20 different geometries 3 different physical quantity in 3 dimensions vectors will be matrix multiplied for a geometry and physics processing neural network layer which will be calculated with genetic selective algorithm which has neural connection coefficient vectors of 5 different groups of population each having 100 different gene sequence (population) and 20 neural nodes. In a such case you may benefit calculating it at once or you can somehow serialize this calculation to two flat arrays and send it to your gpu or cpu according to your concurrent free ram amount. In such case you may want to understand how this calculations work.

In any case of multidimensional matrix calculation combinations of vectors are calculated term by term You may want to multiply first term with seconds vectors last term and sum the rest. It is up to you but understanding how it works is important. So here is a simple illustration I have used to undestand this.

[ [ [A,B],[C,D],[E,F] ], [ [G,H],[I,J],[K,L] ] ] --> ABCDEFGHIJKL [ [ [1,2],[3,4],[5,6] ], [ [7,8],[9,0],[€,$] ] ] --> 1234567890€$

use operator(multiply) term by term shifting first array by amount of vector size (2)

ABCDEFGHIJKL   CDEFGHIJKLAB  FGHIJKLABCDE ...
1234567890€$   1234567890€$  1234567890€$ ...

Here comes all combinations

append all of them and reshape and use another operator (+)

[A+2*B],[C*3+D*4],[E*5,F*6] ...... [I*€+J*$] 

Hope this helps and saves time to grasp this.

Answer 2

If you have numpy arrays you can use the np.dot function for this:

np.dot(A, B)

It will do exactly what you want, i.e. "contract" the last axis of A with the first axis of B:

For 2-D arrays it is equivalent to matrix multiplication, and for 1-D arrays to inner product of vectors (without complex conjugation). For N dimensions it is a sum product over the last axis of a and the second-to-last of b:

 dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])