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.

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]]])

It seems you are doing matrix-multiplications for each slice along the first axis. For the same, you can use np.einsum like so -

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

We can also use np.matmul -

np.matmul(A,B)

On Python 3.x, this matmul operation simplifies with @ operator -

A @ B

Benchmarking

Approaches -

def einsum_based(A,B):
    return np.einsum('ijk,ikl->ijl',A,B)

def matmul_based(A,B):
    return np.matmul(A,B)

def forloop(A,B):
    N = A.shape[0]
    C = np.zeros((N,2,2))
    for i in range(N):
        C[i] = np.dot(A[i], B[i])
    return C

Timings -

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

In [45]: %timeit einsum_based(A,B)
    ...: %timeit matmul_based(A,B)
    ...: %timeit forloop(A,B)
100 loops, best of 3: 3.08 ms per loop
100 loops, best of 3: 3.04 ms per loop
100 loops, best of 3: 10.9 ms per loop

Let's assume b=np.ones((3,2,2)) and a=np.array([1,2,3]). I really do like the answer of @Alok which uses the simple a[:, None, None] * b which surely works with your problem. What I dislike with this formulation is that it's quite dimension specific. What I mean is that it can only be used with 3 dimensional arrays, which was not true in my problem, where b could be a 1D or a 3D array with the exact same length for axis 0. I hence found a way to accommodate it to my problem :

broad_a = np.broadcast_to(a, b.T.shape).T
result = broad_a * b
print(result)
[[
[1,1]
[1,1]]
[
[2,2]
[2,2]]
[
[3,3]
[3,3]]]

Giving also the intended result for your case.

You could do this using np.einsum:

In [142]: old = orig(a,b)

In [143]: new = np.einsum('ijk,ikl->ijl', a, b)

In [144]: np.allclose(old, new)
Out[144]: True

One advantage of using einsum is that you can almost read off what it's doing from the indices: leave the first axis alone (i), and perform a matrix multiplication on the last two (jk,kl->jl)).

The easiest is:

C = numpy.sum(A * B, -1)

I think this might work too:

C = numpy.einsum("...i,...i->...", A, B)