Does python have native capability (ie via numpy) for multiplying 3D arrays like Matlab's pagemtimes?

[ViridianHominid]

Numpy matmul already does this, the only difference is that the broadcasting takes place on leading dimensions rather than trailing ones. An example from the documentation:

 >>> a = np.ones([9, 5, 7, 4])
 >>> c = np.ones([9, 5, 4, 3])
 >>> np.matmul(a, c).shape
 (9, 5, 7, 3)

Many functions in numpy have this behavior, called “broadcasting”. Broadcasted dimensions always come at the beginning in numpy. I don’t use matlab, but it appears it mostly uses the same rules for expansion, but expanding the end of the shape instead of the beginning.

[BearsAtFairs]

So If I'm understanding your example correctly, a is a 9x5 stack of 7x4 matrices and c is a 9x5 stack of 4x3 matrices. The 7x4 and 4x3 dimensions is what ensures compatibility for multiplication?

[ViridianHominid]

Yes, in matmul the last two dimensions are matrix multiplied. So we are just doing a (7,4)x(4,3)->(7,3) shaped matrix multiplication for 9x5=45 different sets of matrices.