From dc7f119b995d9530188280d2b6fe2d55553ff3df Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beat=20K=C3=BCng?= <beat-kueng@gmx.net>
Date: Sat, 21 Aug 2021 19:06:49 +0200
Subject: [PATCH] geninv(): improve runtime performance and reduce stack usage

- use associativity of matrix operations to reduce size of temporary
  matrices and number of multiplications in the M <= N case
- minimize the number of temporary matrices required
---
 matrix/PseudoInverse.hpp | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

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<Type, N, M> geninv(const Matrix<Type, M, N> & G)
         SquareMatrix<Type, M> A = G * G.transpose();
         SquareMatrix<Type, M> L = fullRankCholesky(A, rank);
 
-        SquareMatrix<Type, M> LtL = L.transpose() * L;
+        A = L.transpose() * L;
         SquareMatrix<Type, M> X;
-        if (!inv(LtL, X, rank)) {
+        if (!inv(A, X, rank)) {
             return Matrix<Type, N, M>(); // 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<Type, N> A = G.transpose() * G;
         SquareMatrix<Type, N> L = fullRankCholesky(A, rank);
 
-        SquareMatrix<Type, N> LtL = L.transpose() * L;
+        A = L.transpose() * L;
         SquareMatrix<Type, N> X;
-        if(!inv(LtL, X, rank)) {
+        if(!inv(A, X, rank)) {
             return Matrix<Type, N, M>(); // 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();
     }
 }