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.

As hpaulj suggests in the comments, the easiest way to do this is np.einsum with a dense matrix:

>>> a = np.random.randn(8, 9)
>>> b = np.random.randn(9, 12, 17)
>>> c = np.einsum('ij,jkl->ikl', a, b)
>>> c.shape
(8, 12, 17)

In numpy @ does matrix multiplication

While * does element wise multiplication or Hadamard product

Also, using sum and broadcasting:

np.sum(A * B[:, np.newaxis, :], axis=2)

Einsum seems to be the fastest option though. This would be the slowest since it instantiate the product of elements.