added rpinv (right pseudoinverse)

miscmaths.cc

@@ -519,6 +519,23 @@ namespace MISCMATHS {
       return pinv;
     }
 
+  ReturnMatrix rpinv(const Matrix& mat) // note that mat must be a fat matrix (e.g. x.t() )
+    {
+      // calculates the right psuedo-inverse using SVD
+      Tracer tr("rpinv");
+      DiagonalMatrix D;
+      Matrix U, V;
+      SVD(mat.t(),D,U,V); // feed transpose of input into SVD
+      float tol;
+      tol = MaximumAbsoluteValue(D) * Max(mat.Nrows(),mat.Ncols()) * 1e-16;
+      for (int n=1; n<=D.Nrows(); n++) {
+	if (fabs(D(n,n))>tol)  D(n,n) = 1.0/D(n,n);
+      }
+      Matrix rpinv = U * D * V.t();
+      rpinv.Release();
+      return rpinv;
+    }
+
   int rank(const Matrix& X)
     {
       // calculates the rank of matrix X

miscmaths.h

@@ -120,6 +120,7 @@ namespace MISCMATHS {
   int diag(DiagonalMatrix& m, const ColumnVector& diagvals);
   ReturnMatrix diag(const Matrix& mat);
   ReturnMatrix pinv(const Matrix& mat);
+  ReturnMatrix rpinv(const Matrix& mat);
   int rank(const Matrix& X);
   ReturnMatrix sqrtaff(const Matrix& mat);