diff --git a/matrix/PseudoInverse.hpp b/matrix/PseudoInverse.hpp index 2463529640..a2df21b5d2 100644 --- a/matrix/PseudoInverse.hpp +++ b/matrix/PseudoInverse.hpp @@ -28,23 +28,27 @@ Matrix geninv(const Matrix & G) SquareMatrix A = G * G.transpose(); SquareMatrix L = fullRankCholesky(A, rank); - SquareMatrix LtL = L.transpose() * L; + A = L.transpose() * L; SquareMatrix X; - if (!inv(LtL, X, rank)) { + if (!inv(A, X, rank)) { return Matrix(); // LCOV_EXCL_LINE -- this can only be hit from numerical issues } - return G.transpose() * L * X * X * L.transpose(); + // doing an intermediate assignment reduces stack usage + A = X * X * L.transpose(); + return G.transpose() * (L * A); } else { SquareMatrix A = G.transpose() * G; SquareMatrix L = fullRankCholesky(A, rank); - SquareMatrix LtL = L.transpose() * L; + A = L.transpose() * L; SquareMatrix X; - if(!inv(LtL, X, rank)) { + if(!inv(A, X, rank)) { return Matrix(); // LCOV_EXCL_LINE -- this can only be hit from numerical issues } - return L * X * X * L.transpose() * G.transpose(); + // doing an intermediate assignment reduces stack usage + A = X * X * L.transpose(); + return (L * A) * G.transpose(); } }