2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Is it possible to elaborate a bit more about what do you mean "to increase numerical stability"? $\endgroup$ – Alone Programmer Dec 10 '19 at 17:44
  • 1
    $\begingroup$ The latter option. $\endgroup$ – user1767774 Dec 11 '19 at 7:27
  • 1
    $\begingroup$ In that case, your matrices aren't orthogonal matrices and the answer by @whpowell96 isn't relevant. $\endgroup$ – Brian Borchers Dec 11 '19 at 13:58
  • 1
    $\begingroup$ You're starting with a 300 dimensional space, and each projection would reduce this dimension by at least one. Unless there's significant repition in this sequence of projections, you'd be left with only the 0 vector after 400 projections. Could you explain where this problem is coming from? $\endgroup$ – Brian Borchers Dec 11 '19 at 16:08
  • 1
    $\begingroup$ Yes, in practice I work with matrices of different dimensions, and the number of products changes accordingly; 300x300 matrices are multiplied about 100 times. The problem comes from a research in ML, and each projection matrix in the product should remove a different kind of information from the input vectors. $\endgroup$ – user1767774 Dec 11 '19 at 19:29
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks. Is there also something to do with regard to general rounding errors - like working in log space as is commonly done when multiplying probabilities? $\endgroup$ – user1767774 Dec 10 '19 at 19:53
  • 1
    $\begingroup$ The loss of orthogonality is a result of accumulated round off errors. This should be the only significant source of error since orthogonal matrices are well-conditioned $\endgroup$ – whpowell96 Dec 10 '19 at 20:25
  • $\begingroup$ Thanks. Are there general good practices for the general case, in which the matrices are not orthogonal? is it possible to perform the product in log space? $\endgroup$ – user1767774 Dec 10 '19 at 20:52
  • 2
    $\begingroup$ Not that I am aware of. In the general case, the errors get scaled further by the condition numbers of the factors. The multiplication may be more stable if each factor is stored in some decomposition like QR or SVD but that would be inefficient to compute for this many factors $\endgroup$ – whpowell96 Dec 11 '19 at 3:58
  • 2
    $\begingroup$ The OP's matrices are not "orthogonal matrices" Rather they are the matrices of orthogonal projection operations. $\endgroup$ – Brian Borchers Dec 11 '19 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.