I have to calculate in numpy the matrix-product of many matrices (~400). Are there common practices to increase numerical stability?
If this is relevant, the matrices are $300\times 300$ orthogonal projection matrices.
Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal part. Since your matrices are square, this will cost $O(n^3)$, so comparable to the matrix multiplications.