Not exactly sure, what you are trying to achieve. Maybe you could give an example of your input and your expected output. One possibility is:

import numpy as np

x = np.array([[1, 2], [1, 2], [1, 2]])
y = np.array([1, 2, 3])
res = x * np.transpose(np.array([y,]*2))

This will multiply each column of x with y, so the result of the above example is:

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

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.

Note that np.dot behaves differently with different dimension arrays, from here:

For 2-D arrays it is equivalent to matrix multiplication,

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

In your example z has 3 dimensions.

Also, there is an important distinction between a numpy array and a numpy matrix:

You are using arrays in your example. Numpy matrices are 2D only and matrix multiplication is achived with the * operator. (With arrays this perfroms element by element multiplication)

Ha, it can be done in just one line: np.einsum('nmk,nkj->mj',A,B).

See Einstein summation: http://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html

Not the same problem but the idea is quite much the same, see discussions and alternative methods in this topic we just discussed: numpy multiply matrices preserve third axis

Don't name your variable sum, you override the build-in sum.

As @Jaime pointed out, the loop is actually faster for dimensions of these size. In fact a solution based on map and sum is, albeit simpler, even slower:

In [19]:

%%timeit
SUM = np.zeros([20,5])
for i in range(len(A)):
  SUM += np.dot(A[i],B[i])
10000 loops, best of 3: 115 µs per loop
In [20]:

%timeit np.array(map(np.dot, A,B)).sum(0)
1000 loops, best of 3: 445 µs per loop
In [21]:

%timeit np.einsum('nmk,nkj->mj',A,B)
1000 loops, best of 3: 259 µs per loop

Thing are different with larger dimension:

n_examples = 1000
A = np.random.randn(n_examples, 20,1000)
B = np.random.randn(n_examples, 1000,5)

And:

In [46]:

%%timeit
SUM = np.zeros([20,5])
for i in range(len(A)):
  SUM += np.dot(A[i],B[i])
1 loops, best of 3: 191 ms per loop
In [47]:

%timeit np.array(map(np.dot, A,B)).sum(0)
1 loops, best of 3: 164 ms per loop
In [48]:

%timeit np.einsum('nmk,nkj->mj',A,B)
1 loops, best of 3: 451 ms per loop