From 38884acf3a315f63e2e5f1e9948a5548a0560cab Mon Sep 17 00:00:00 2001
From: Stephen Smith <steve@fmrib.ox.ac.uk>
Date: Mon, 13 Feb 2006 11:32:18 +0000
Subject: [PATCH] added rpinv (right pseudoinverse)
---
miscmaths.cc | 17 +++++++++++++++++
miscmaths.h | 1 +
2 files changed, 18 insertions(+)
diff --git a/miscmaths.cc b/miscmaths.cc
index 83270df..7992bf9 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 4452b44..6474e96 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);
--
GitLab