Matrix multiplication over specific dimensions in tensorflow (or numpy)

Question

I hope I'm not missing anything obvious here, but I've scoured the inter-webs to no avail, and finally come to ask here...
Here's a really dry and simple description of what I'd like to do:
Say I've got a tensor of shape (20, 40, 3, 5), and another tensor of shape (20, 40, 5, 7). The first two dimension sizes are to be kept as are, and are purposely identical for the two tensors. The last two dimensions on the other hand, are to be (matrix-)multiplied, matmul style. Meaning my resulting tensor would be of shape (20, 40, 3, 7). How can that be done??
I realize I can theoretically just loop over the first two dimensions and use tf.matmul() directly, but that's an absolute no-go due to runtime, efficiency, model-trainer and GPU world-wide protests, and my conscience if that's of any weight :-).

I've unfortunately disregarded as "not what I need" the following options:
tf.tensordot would give me an output of shape (20, 40, 3, 20, 40, 7). No good.
tf.scan is only good for the first dimension if I'm reading it correctly (suitable for RNNs maybe? Not my case anyhow).
tf.matmul works for tensors of rank >= 2, but works like @ over the last and first dimensions respectively. Again, not my case.

So again - how can this be done?
A numpy answer that helps me get in the right direction would also be very helpful, but I'll need a tf implementation at the end of the day.

Thanks in advance, and sorry if I'm missing something dumb.

---

The following is closer to what I need, but less clear and so is being written separately:
The first two dimensions are spatial dimensions of an image. The last two are actually square matrices, obtained via tf.contrib.distributions.fill_triangular, and are being multiplied (along with an appropriate transpose on one of them) to obtain covariance matrices associated to each spatial coordinate. I don't know if that helps in any way, but it gives some context at the very least. Also, there might or might not be a batch dimension as well, but I'm assuming that solving the 4-D tensor case would be generalizable enough.

Answer 1

Posting this for future reference:
From numpy matmul doc:

If either argument is N-D, N > 2, it is treated as a stack of matrices residing in the last two indexes and broadcast accordingly.

For dimensions >2 it will treat it as a stack of matrices, attempting to matmul the last 2 dimensions, resulting with a np array as the OP required.
For example:

import numpy as np

A = np.ones((1,2,1,2))
B = np.ones((1,2,2,1))

print(A.shape)
print(B.shape)
print(np.matmul(A,B).shape)

with result:

(1, 2, 1, 2)
(1, 2, 2, 1)
(1, 2, 1, 1)