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)

We can use np.einsum -

# a, b are 3D and 2D arrays respectively
np.einsum('ijk,lk->ilj', a, b)

Alternatively, with np.matmul/@-operator on Python3.x -

np.matmul(a,b.T[None]).swapaxes(1,2)

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,

Yes that's possible. If you have mutiple batch dimensions in both operatns, you can use the broadcasting. In this case the last two dimensions of each operand are interpreted as a matrix size. (I recommend looking it up in the documentation.)

So you need an additional dimension for your vectors b, to make them a n x 1 "matrix" (column vector):

# original implementation
b1 = b.unsqueeze(1)
r1 = z * b1
r1 = torch.sum(r1, dim=-1)
print(r1.shape)

# using torch.matmul
r2 = torch.matmul(z, b.unsqueeze(2))[...,0]
print(r2.shape)
print((r1-r2).abs().sum())  # should be zero if we do the same operation

Alternatively, torch.einsum also makes this very straightforward.

# using torch.einsum
r3 = torch.einsum('ijk,ik->ij', z, b)
print((r1-r3).abs().sum())  # should be zero if we do the same operation

einsum is a very powerful operation that can do a lot of things: you can permute tensor dimensions, sum along them, or perform scalar products, all with or without broadcasting. It is derived from the Einstein summation convention mostly used in physics. The rough idea is that you give every dimension of your operans a name, and then, using these names define what the output should look like. I think it is best to read the documentation. In our case we have a 4 x 3 x 2 tensor as well as a 4 x 2 tensor. So the let's call the dimensions of the first tensor ijk. Here i and k should be considered the same as the dimensions of the second tensor, so this one can be described as ik. Finally the output should have clearly be ij (it mus be a 4 x 3 tensor). From this "signature" ijk, ik -> ij it is clear that the dimension i is preserved, and the dimensions k must be "summe/multiplied" away (scalar product).

You can achieve that like this.

L --> [?, M, M]
T --> [M, N]
tensordot(L,T) axes [[2], [0]] --> [?,M, N]
F --> [?, N] --> expand axis --> [?, N, 1]
matmul [?, M, N], [?, N, 1] --> [?, M, 1] --> squeeze --> [?, M]

putting together:

tf.squeeze(tf.matmul(tf.tensordot(L,T, axes=[[2],[0]]),F[...,None]))

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)

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.