I discovered np.einsum, the most useful of all Numpy library functions.
Suppose you want to compute:
$$ T^{svi} = \sum_x \sum_j P^{sx} M^{v}_{xj} Q^{ij} $$
Then you can just describe this operation using indices:
T = np.einsum("sx,vxj,ij -> svi", P, M, Q)
It's more compact than the LaTeX expression!
Another example:
$$
M^{v}_{xj} = \sum_{x} \sum_{j} T^{svi} P^{-1}_{sx} Q^{-1}_{ij}
$$
becomes:
M = np.einsum("svi,sx,ij -> vxj", T, P_inv, Q_inv)