diff --git a/miscmaths.cc b/miscmaths.cc
index 83270df14515fb7022471060c498b3c85affcc6d..7992bf9d0e56f001f06da2015406919857fee8a2 100644
--- a/miscmaths.cc
+++ b/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
diff --git a/miscmaths.h b/miscmaths.h
index 4452b44c34ae33b305c3da23913be06c14a511c6..6474e963bf12b864133339aef65962e4b6d290ea 100644
--- a/miscmaths.h
+++ b/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);