It seems you are trying to lose the last two axes of the first array against the only two axes of the second weight array with that matrix-multiplication. We could translate that idea into NumPy code with np.tensordot and assuming arr1 and arr2 as the input arrays respectively, like so -

np.tensordot(arr1,arr2,axes=([1,2],[0,1]))

Another simpler way to put into NumPy code would be with np.einsum, like so -

np.einsum('ijk,jk',arr1,arr2)

Use np.tensordot and then swap axes. So, use one of these -

np.tensordot(B,A,axes=((1),(1))).swapaxes(0,1)
np.tensordot(A,B,axes=((1),(1))).swapaxes(1,2)

We can reshape A to 2D after swapping axes, use 2D matrix multiplication with np.dot and reshape and swap axes to maybe gain marginal performance boost.

Timings -

# Original approach
def orgapp(A,B):
    m = A.shape[0]
    n = B.shape[0]
    r = A.shape[2]
    c = np.zeros((m,n,r))
    for i in range(len(A)):
        c[i] = np.dot(B, A[i,:,:])
    return c  

In [91]: n = 10000
    ...: A = np.random.rand(n,2,2)
    ...: B = np.random.rand(2,2)

In [92]: %timeit orgapp(A,B)
100 loops, best of 3: 12.2 ms per loop

In [93]: %timeit np.tensordot(B,A,axes=((1),(1))).swapaxes(0,1)
1000 loops, best of 3: 191 µs per loop

In [94]: %timeit np.tensordot(A,B,axes=((1),(1))).swapaxes(1,2)
1000 loops, best of 3: 208 µs per loop

# @Bitwise's solution
In [95]: %timeit np.flip(np.dot(A,B).transpose((0,2,1)),1)
1000 loops, best of 3: 697 µs per loop

You can use np.tensordot -

np.tensordot(A,B,axes=((2),(1)))

Related post to understand tensordot.

einsum equivalent would be -

np.einsum('ijk,lk->ijl', A, B)

We can also use A.dot(B.T), but that would be looping under the hoods. So, might not be the most preferred one, but it's a compact solution,

The problem is that numpy views multidimensional arrays as stacks of matrices, and always the last two dimensions are assumed to be the linear space dimensions. This means that the dot product will not work by collapsing the first dimension of your 3d array.

Instead the simplest thing you can do is to reshape your 3d array into a 2d one, doing the matrix multiplication, and reshaping back into a 3d array. This will also make use of optimised BLAS code which is one of the great advantages of numpy.

import numpy as np 

S_pinv = np.random.rand(3, 4)
images = np.random.rand(4, 5, 6)

# error: 
# (S_pinv @ images).shape 
res_shape = S_pinv.shape[:1] + images.shape[1:]  # (3, 5, 6) 
res = (S_pinv @ images.reshape(images.shape[0], -1)).reshape(res_shape)
print(res.shape)  # (3, 5, 6)

So instead of (3,n) x (n,h,w) we do (3,n) x (n, h*w) -> (3, h*w) which we reshape back to (3, h, w). Reshaping is free, because this doesn't mean any actual manipulation of data in memory (only a reinterpretation of the single block of memory that underlies the array), and as I said proper matrix products are highly optimized in numpy.

Since you asked for a more intuitive way, here's an alternative making use of numpy.einsum. It will probably be slower, but it's very transparent if you get a little bit used to its notation:

res_einsum = np.einsum('tn,nhw -> thw', S_pinv, images)
print(np.array_equal(res, res_einsum))  # True

This notation names each of the dimensions of the input arrays: for S_pinv the first and second dimensions are named t and n, respectively, and similarly n, h and w for images. The output is set to have dimensions thw which means that any remaining dimensions that are not present in the output shape will be summed along after multiplying the input arrays. This is exactly what you need.

As you noted in a comment, you could also transpose your arrays so that np.dot finds the right dimensions in the right place. But this will also be slow because this might lead to copies in memory, or at least suboptimal looping over your arrays.

I made a quick timing comparison using the following defininitions:

def reshaped(S_pinv, images): 
    res_shape = S_pinv.shape[:1] + images.shape[1:] 
    return (S_pinv @ images.reshape(images.shape[0], -1)).reshape(res_shape)

def einsummed(S_pinv, images): 
    return np.einsum('tn,nhw -> thw', S_pinv, images) 

def transposed(S_pinv, images): 
    return (S_pinv @ images.transpose(2, 0, 1)).transpose(1, 2, 0)          

And here's the timing test using IPython's built-in %timeit magic, and some more realistic array sizes:

>>> S_pinv = np.random.rand(3, 100) 
... images = np.random.rand(100, 200, 300) 
... args = S_pinv, images 
... %timeit reshaped(*args) 
... %timeit einsummed(*args) 
... %timeit transposed(*args)                                          
5.92 ms ± 460 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
15.9 ms ± 190 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
44.5 ms ± 329 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

For array of dimension 3 (or, a rank-3 tensor) that you have, you can use np.einsum doc for more complex matrix multiplications. In your particular case, you can use the following

>>> import numpy as np
>>> x = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> y = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> np.einsum('ijk,ikl->ijl', x, y)  # still shape (3, 3, 3)

In particular, the einsum expression 'ijk,ikl->ijl' means for each ith matrix, do a regular matrix multiplication jk,kl->jl and put the result in the ith entry in the resulting tensor (ndarray). A more general form of this process could be

np.einsum('...jk,...kl->...jl', x, y)

where you can have arbitrary number of dimensions in front of each tensor (ndarray) that you have.

See the following for a complete example:

>>> import numpy as np
>>> x = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> x
array([[[0, 0, 1],
        [2, 2, 1],
        [2, 1, 1]],

       [[2, 0, 2],
        [2, 2, 1],
        [2, 2, 2]],

       [[2, 2, 2],
        [1, 1, 2],
        [0, 2, 2]]])
>>> y = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> y
array([[[0, 0, 1],
        [2, 1, 0],
        [0, 0, 2]],

       [[1, 2, 0],
        [2, 0, 1],
        [2, 2, 1]],

       [[0, 2, 1],
        [0, 1, 0],
        [0, 2, 1]]])
>>> np.einsum('ijk,ikl->ijl', x, y)
array([[[ 0,  0,  2],
        [ 4,  2,  4],
        [ 2,  1,  4]],

       [[ 6,  8,  2],
        [ 8,  6,  3],
        [10,  8,  4]],

       [[ 0, 10,  4],
        [ 0,  7,  3],
        [ 0,  6,  2]]])
>>> np.einsum('...ij,...jk->...ik', x, y)
array([[[ 0,  0,  2],
        [ 4,  2,  4],
        [ 2,  1,  4]],

       [[ 6,  8,  2],
        [ 8,  6,  3],
        [10,  8,  4]],

       [[ 0, 10,  4],
        [ 0,  7,  3],
        [ 0,  6,  2]]])

Sorry for the necromancy, but this answer can be substantially improved upon, using the invaluable np.einsum.

import numpy as np

D,M,N,R = 1,2,3,4
A = np.random.rand(M,N,R)
B = np.random.rand(N,D,R)

print np.einsum('mnr,ndr->mdr', A, B).shape

Note that it has several advantages: first of all, its fast. np.einsum is well-optimized generally, but moreover, np.einsum is smart enough to avoid the creation of an MxNxR temporary array, but performs the contraction over N directly.

But perhaps more importantly, its very readable. There is no doubt that this code is correct; and you could make it a lot more complicated without any trouble.

Note that the dummy 'D' axis can simply be dropped from B and the einsum statement if you wish.