pseudoinverse

Former-commit-id: ce0bc49c10ab9cfa1be71ea1bf7541f4d2f70780

include/igl/pinv.cpp

@@ -0,0 +1,35 @@
+#include "pinv.h"
+#include <limits>
+#include <cmath>
+
+template <typename DerivedA, typename DerivedX>
+void igl::pinv(
+  const Eigen::MatrixBase<DerivedA> & A,
+  typename DerivedA::Scalar tol,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  Eigen::JacobiSVD<DerivedA> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV );
+  typedef typename DerivedA::Scalar Scalar;
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & U = svd.matrixU();
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> & V = svd.matrixV();
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,1> & S = svd.singularValues();
+  if(tol < 0)
+  {
+    const Scalar smax = S.array().abs().maxCoeff();
+    tol = 
+      (Scalar)(std::max(A.rows(),A.cols())) *
+      (smax-std::nextafter(smax,std::numeric_limits<Scalar>::epsilon()));
+  }
+  const int rank = (S.array()>0).count();
+  X = (V.leftCols(rank).array().rowwise() * 
+      (1.0/S.head(rank).array()).transpose()).matrix()*
+    U.leftCols(rank).transpose();
+}
+
+template <typename DerivedA, typename DerivedX>
+void igl::pinv(
+  const Eigen::MatrixBase<DerivedA> & A,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  return pinv(A,-1,X);
+}

include/igl/pinv.h

@@ -0,0 +1,29 @@
+#ifndef IGL_PINV_H
+#define IGL_PINV_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute the Moore-Penrose pseudoinverse
+  //
+  // Inputs:
+  //   A  m by n matrix 
+  //   tol  tolerance (if negative then default is used)
+  // Outputs:
+  //   X  n by m matrix so that A*X*A = A and X*A*X = X and A*X = (A*X)' and
+  //     (X*A) = (X*A)'
+  template <typename DerivedA, typename DerivedX>
+  void pinv(
+    const Eigen::MatrixBase<DerivedA> & A,
+    typename DerivedA::Scalar tol,
+    Eigen::PlainObjectBase<DerivedX> & X);
+  // Wrapper using default tol
+  template <typename DerivedA, typename DerivedX>
+  void pinv(
+    const Eigen::MatrixBase<DerivedA> & A,
+    Eigen::PlainObjectBase<DerivedX> & X);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "pinv.cpp"
+#endif
+#endif