New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 257 - pseudoinverse
pseudoinverse
Status: CONFIRMED
Product: Eigen
Classification: Unclassified
Component: SVD
3.0
All All
: Low Feature Request
Assigned To: Nobody
:
: 412 (view as bug list)
Depends on: ExprEval
Blocks:
  Show dependency treegraph
 
Reported: 2011-04-28 15:52 UTC by Cyrille Berger
Modified: 2016-07-04 14:40 UTC (History)
12 users (show)



Attachments

Description Cyrille Berger 2011-04-28 15:52:16 UTC
It would be nice to have a build in pseudoinverse.

In the meantime, some snipets are available:

http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2010/01/msg00187.html
Comment 1 Cyrille Berger 2011-04-28 16:06:35 UTC
Here is an updated snipet

#include <Eigen/SVD>

template<typename _Matrix_Type_>
bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
{
  if(a.rows()<a.cols())
      return false;

  Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd();

  typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs().maxCoeff();

  result = svd.matrixV() * _Matrix_Type_(_Matrix_Type_( (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().
      array().inverse(), 0) ).diagonal()) * svd.matrixU().adjoint();
}


(I had two add two temporary matrixes, since it seems select does not have the diagonal function and multiplying the matrix with the diagonal matrix object was triggering a compile error... but I don't have the time right now to investigate more those issues)
Comment 2 Jitse Niesen 2011-04-28 16:32:28 UTC
Parts of this can ideally be used to resolve bug 241 (using LDLT to solve linear equations with singular matrices).
Comment 3 Cyrille Berger 2011-04-28 16:37:09 UTC
Bah too fast... Here is a version tested when executing, with one less temporary:

template<typename _Matrix_Type_>
bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
{
  if(a.rows()<a.cols())
      return false;

  Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);

  typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs().maxCoeff();
  
  result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().
      array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
}
Comment 4 Christoph Hertzberg 2011-04-28 18:49:51 UTC
(In reply to comment #3)
> Bah too fast... Here is a version tested when executing, with one less
> temporary:
> 
> template<typename _Matrix_Type_>
> bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double
> epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
> {
>   if(a.rows()<a.cols())
>       return false;

Why this restriction? (If this is due to a JacobiSVD restriction, one can use pinv(A) == pinv(A.transpose()).transpose()). That way pseudoInverse would always succeed.

>   Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
> Eigen::ComputeFullV);

I think you just need thin U and V.


>   typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
> a.rows()) * svd.singularValues().array().abs().maxCoeff();

I assume singularValues are sorted and non-negative, so the last factor could be replaced by svd.singularValues()[0].

>   result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() >
> tolerance).select(svd.singularValues().
>       array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
> }


It would be great, if JacobiSVD::solve had an additional parameter epsilon; that way pseudoInverse(a) would be just:

a.jacobiSvd(ComputeThinU | ComputeThinV).solve(MatrixXd::Identity(a.rows()), epsilon);
Comment 5 Philippe Hamelin 2011-05-19 15:39:35 UTC
My own pseudoinverse looks like this and yes, this code will work only with ThinU and ThinV.
Comment 6 catriel.carbonera 2012-01-30 02:07:49 UTC
*** Bug 412 has been marked as a duplicate of this bug. ***
Comment 7 hsieh.boss 2012-04-24 12:13:02 UTC
(In reply to comment #3)
> Bah too fast... Here is a version tested when executing, with one less
> temporary:
> 
> template<typename _Matrix_Type_>
> bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double
> epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
> {
>   if(a.rows()<a.cols())
>       return false;
> 
>   Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
> Eigen::ComputeFullV);
> 
>   typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
> a.rows()) * svd.singularValues().array().abs().maxCoeff();
> 
>   result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() >
> tolerance).select(svd.singularValues().
>       array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
> }

I try to use this function, but it will cause assert failed when a is not square matrix.
I have to create a non square matrix to save the sigma, like

_Matrix_Type_ sigma = _Matrix_Type_::Zero(a.cols(), a.rows());
sigma.block(0, 0, a.cols(), a.cols())<<_Matrix_Type_(_Matrix_Type_( (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0) ).asDiagonal());

result = (svd.matrixV() * sigma) * svd.matrixU().adjoint();

Is there any better method to avoid this problem?
I am just start to use Eigen recently.
Comment 8 Koubi 2012-05-09 22:28:49 UTC
(In reply to comment #7)
> (In reply to comment #3)
> > Bah too fast... Here is a version tested when executing, with one less
> > temporary:
> > 
> > template<typename _Matrix_Type_>
> > bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double
> > epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
> > {
> >   if(a.rows()<a.cols())
> >       return false;
> > 
> >   Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
> > Eigen::ComputeFullV);
> > 
> >   typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
> > a.rows()) * svd.singularValues().array().abs().maxCoeff();
> > 
> >   result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() >
> > tolerance).select(svd.singularValues().
> >       array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
> > }
> 
> I try to use this function, but it will cause assert failed when a is not
> square matrix.
> I have to create a non square matrix to save the sigma, like
> 
> _Matrix_Type_ sigma = _Matrix_Type_::Zero(a.cols(), a.rows());
> sigma.block(0, 0, a.cols(), a.cols())<<_Matrix_Type_(_Matrix_Type_(
> (svd.singularValues().array().abs() >
> tolerance).select(svd.singularValues().array().inverse(), 0) ).asDiagonal());
> 
> result = (svd.matrixV() * sigma) * svd.matrixU().adjoint();
> 
> Is there any better method to avoid this problem?
> I am just start to use Eigen recently.

You can just change "Eigen::ComputeFullU | Eigen::ComputeFullV" in "Eigen::ComputeThinU |Eigen::ComputeThinV". I tried it and it works.
Comment 9 Rafael 2012-06-19 20:04:05 UTC
Hi everyone,

I started today using Eigen and my purpose is compute a pseudoinverse of 80,000 x 80,000. I wrote the code below and I would like to know if someone can help me to solve this problem. Basically, I don't know how to build the loop to receive the function suggested below.:

Thanks in advance,

Rafael

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string>
#include <string.h>
#include <vector>
#include <sstream>
#include <algorithm>
#include <cstdlib>
#include <ctime>


#include <sys/types.h>
#include <unistd.h>

#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/SVD>

typedef Eigen::Matrix< float_t, Eigen::Dynamic, Eigen::Dynamic > MatrixXf;


using namespace std;
using namespace Eigen;
using Eigen::MatrixXd;
using Eigen::MatrixXf;
using Eigen::Matrix3d;



int main()
{
	ifstream MA;
	MA.open("A.out");

	int dimension = 8;

	// TEST EIGEN
	MatrixXf m;
	m = MatrixXf::Zero(dimension, dimension);

	if (MA.is_open())
	{
		while(!MA.eof())
		{
			for(int i=0;i<dimension;i++)
				for(int j=0;j<dimension;j++)
					MA >> m(i,j);
		}
		MA.close();
	}
	else
		cout << "Unable to open file IN_AOUT";

	std::cout << "OK!" << std::endl;

	Eigen::JacobiSVD <MatrixXf>svd(m,ComputeThinU | ComputeThinV);


	cout << "Its singular values are:" << endl <<  svd.singularValues() << endl; // S
	cout << "Its left singular vectors are the columns of the thin U matrix:" << endl << svd.matrixU() << endl; // U
	cout << "Its right singular vectors are the columns of the thin V matrix:" << endl << svd.matrixV() << endl; // V

	cout << svd.matrixV()*svd.singularValues().head().().inverse()*svd.matrixU().transpose();


    return 0;
}








(In reply to comment #8)
> (In reply to comment #7)
> > (In reply to comment #3)
> > > Bah too fast... Here is a version tested when executing, with one less
> > > temporary:
> > > 
> > > template<typename _Matrix_Type_>
> > > bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double
> > > epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
> > > {
> > >   if(a.rows()<a.cols())
> > >       return false;
> > > 
> > >   Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
> > > Eigen::ComputeFullV);
> > > 
> > >   typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
> > > a.rows()) * svd.singularValues().array().abs().maxCoeff();
> > > 
> > >   result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() >
> > > tolerance).select(svd.singularValues().
> > >       array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
> > > }
> > 
> > I try to use this function, but it will cause assert failed when a is not
> > square matrix.
> > I have to create a non square matrix to save the sigma, like
> > 
> > _Matrix_Type_ sigma = _Matrix_Type_::Zero(a.cols(), a.rows());
> > sigma.block(0, 0, a.cols(), a.cols())<<_Matrix_Type_(_Matrix_Type_(
> > (svd.singularValues().array().abs() >
> > tolerance).select(svd.singularValues().array().inverse(), 0) ).asDiagonal());
> > 
> > result = (svd.matrixV() * sigma) * svd.matrixU().adjoint();
> > 
> > Is there any better method to avoid this problem?
> > I am just start to use Eigen recently.
> 
> You can just change "Eigen::ComputeFullU | Eigen::ComputeFullV" in
> "Eigen::ComputeThinU |Eigen::ComputeThinV". I tried it and it works.
Comment 10 Christoph Hertzberg 2013-08-23 16:19:07 UTC
Once this is solved, this FAQ entry can be changed:
http://eigen.tuxfamily.org/index.php?title=FAQ#Is_there_a_method_to_compute_the_.28Moore-Penrose.29_pseudo_inverse_.3F
Comment 11 Ojaswa 2013-10-28 12:57:50 UTC
I would like to compute pseudoinverse of a complex matrix. My function looks like:

template<typename _Matrix_Type_>
bool pseudoInverse(const _Matrix_Type_ &a, _Matrix_Type_ &result, double epsilon = std::numeric_limits<double>::epsilon())
{
  if(a.rows() < a.cols())
   return false;
  Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV);
  double tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs().maxCoeff();
  result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0) ).asDiagonal() * svd.matrixU().adjoint();
}

When I apply this function to matrix type MatrixXcd, I get a lot of compile errors (perhaps in last line 'result = ...'). Please see errors below. How can I modify it to work for complex matrices?

Many Thanks.
-Ojaswa
------------------------------------------ERRORS-------------------------------------------------
$ g++ -I/cygdrive/d/Eigen-3.2/ pseudo_inverse.cpp -o pseudo_inverse
In file included from /cygdrive/d/Eigen-3.2/Eigen/Core:284:0,
                 from /cygdrive/d/Eigen-3.2/Eigen/Dense:1,
                 from pseudo_inverse.cpp:2:
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Assign.h: In member function ‘Derived& Eigen::DenseBase<Derived>::lazyAssign(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, Derived = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’:
/cygdrive/d/Eigen-3.2/Eigen/src/Core/PlainObjectBase.h:411:46:   instantiated from ‘Derived& Eigen::PlainObjectBase<Derived>::lazyAssign(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, Derived = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Assign.h:520:123:   instantiated from ‘static Derived& Eigen::internal::assign_selector<Derived, OtherDerived, false, false>::run(Derived&, const OtherDerived&) [with Derived = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>, OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >]’
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Assign.h:564:89:   instantiated from ‘Derived& Eigen::MatrixBase<Derived>::operator=(const Eigen::DenseBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, Derived = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’
/cygdrive/d/Eigen-3.2/Eigen/src/Core/PlainObjectBase.h:450:7:   instantiated from ‘Derived& Eigen::PlainObjectBase<Derived>::operator=(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, Derived = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Matrix.h:184:35:   instantiated from ‘Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::operator=(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, _Scalar = std::complex<double>, int _Rows = -0x000000001, int _Cols = -0x000000001, int _Options = 0, int _MaxRows = -0x000000001, int _MaxCols = -0x000000001, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Matrix.h:310:7:   instantiated from ‘Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::Matrix(const Eigen::EigenBase<OtherDerived>&) [with OtherDerived = Eigen::Select<Eigen::CwiseUnaryOp<std::binder2nd<std::greater<double> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, -0x000000001, 1> > > > >, _Scalar = std::complex<double>, int _Rows = -0x000000001, int _Cols = -0x000000001, int _Options = 0, int _MaxRows = -0x000000001, int _MaxCols = -0x000000001]’
pseudo_inverse.cpp:15:3:   instantiated from ‘bool pseudoInverse(const _Matrix_Type_&, _Matrix_Type_&, double) [with _Matrix_Type_ = Eigen::Matrix<std::complex<double>, -0x000000001, -0x000000001>]’
pseudo_inverse.cpp:35:35:   instantiated from here
/cygdrive/d/Eigen-3.2/Eigen/src/Core/Assign.h:493:3: error: ‘YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY’ is not a member of ‘Eigen::internal::static_assertion<false>’
------------------------------------------ERRORS-------------------------------------------------
Comment 12 Christoph Hertzberg 2013-10-28 13:24:34 UTC
(In reply to comment #11)
>   result = svd.matrixV() * _Matrix_Type_( (svd.singularValues().array().abs() ...

change that line to:

  result = svd.matrixV() * ( (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0) ).matrix().asDiagonal() * svd.matrixU().adjoint();

The problem is that singularValues() returns an array of real values which will not be converted to a complex matrix automatically. Products between real and complex matrices are implemented, however. Alternatively, you could manually cast it it to complex.
Comment 13 Ojaswa 2013-10-29 09:18:23 UTC
Thanks Christoph for the hint! It works now.
Comment 14 davidhigh 2014-04-08 00:54:02 UTC
I just collected the above codepieces and the suggestions for improvement by Christoph, and here is what I got. Notably,
(i) this code seem to work for any number of rows/columns (no need for the restriction cols<rows)
(ii) it directly returns an Eigen MatrixType instead of changing a reference. Given return value optimization, this should give no overhead.


------------------------------------------------------------
#include<iostream>
#include<Eigen/Core>
#include<Eigen/SVD>


template<typename _Matrix_Type_>
_Matrix_Type_ pseudoInverse(const _Matrix_Type_ &a, double epsilon = std::numeric_limits<double>::epsilon())
{
	Eigen::JacobiSVD< _Matrix_Type_ > svd(a ,Eigen::ComputeThinU | Eigen::ComputeThinV);
	double tolerance = epsilon * std::max(a.cols(), a.rows()) *svd.singularValues().array().abs()(0);
	return svd.matrixV() *  (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint();
}


int main()
{
	Eigen::MatrixXd A(2,3);
	A<< 1, 2, 3, 4, 5, 7;
	std::cout<<A<<std::endl<<std::endl;
	std::cout<<pseudoInverse(A)<<std::endl;

	Eigen::MatrixXd B(3,2);
	B<< 1, 2, 3, 4, 5, 7;
	std::cout<<B<<std::endl<<std::endl;
	std::cout<<pseudoInverse(B)<<std::endl;

	return 0;
}
------------------------------------------------------------


Here is the result of the main() function:

------------------------------------------------------------
1 2 3
4 5 7

     -1.42857           0.6
     0.142857 -2.35922e-016
     0.714286          -0.2

1 2
3 4
5 7

    -2.07143     0.785714     0.142857
         1.5         -0.5 8.46545e-016
------------------------------------------------------------

which agrees with the one I get from this site: http://www.bluebit.gr/matrix-calculator

Best regards,
David
Comment 15 salan668 2014-04-14 09:05:53 UTC
(In reply to comment #14)
> I just collected the above codepieces and the suggestions for improvement by
> Christoph, and here is what I got. Notably,
> (i) this code seem to work for any number of rows/columns (no need for the
> restriction cols<rows)
> (ii) it directly returns an Eigen MatrixType instead of changing a reference.
> Given return value optimization, this should give no overhead.
> 
> 
> ------------------------------------------------------------
> #include<iostream>
> #include<Eigen/Core>
> #include<Eigen/SVD>
> 
> 
> template<typename _Matrix_Type_>
> _Matrix_Type_ pseudoInverse(const _Matrix_Type_ &a, double epsilon =
> std::numeric_limits<double>::epsilon())
> {
>     Eigen::JacobiSVD< _Matrix_Type_ > svd(a ,Eigen::ComputeThinU |
> Eigen::ComputeThinV);
>     double tolerance = epsilon * std::max(a.cols(), a.rows())
> *svd.singularValues().array().abs()(0);
>     return svd.matrixV() *  (svd.singularValues().array().abs() >
> tolerance).select(svd.singularValues().array().inverse(),
> 0).matrix().asDiagonal() * svd.matrixU().adjoint();
> }
> 
> 
> int main()
> {
>     Eigen::MatrixXd A(2,3);
>     A<< 1, 2, 3, 4, 5, 7;
>     std::cout<<A<<std::endl<<std::endl;
>     std::cout<<pseudoInverse(A)<<std::endl;
> 
>     Eigen::MatrixXd B(3,2);
>     B<< 1, 2, 3, 4, 5, 7;
>     std::cout<<B<<std::endl<<std::endl;
>     std::cout<<pseudoInverse(B)<<std::endl;
> 
>     return 0;
> }
> ------------------------------------------------------------
> 
> 
> Here is the result of the main() function:
> 
> ------------------------------------------------------------
> 1 2 3
> 4 5 7
> 
>      -1.42857           0.6
>      0.142857 -2.35922e-016
>      0.714286          -0.2
> 
> 1 2
> 3 4
> 5 7
> 
>     -2.07143     0.785714     0.142857
>          1.5         -0.5 8.46545e-016
> ------------------------------------------------------------
> 
> which agrees with the one I get from this site:
> http://www.bluebit.gr/matrix-calculator
> 
> Best regards,
> David

When I first test these code, it works. But today  it can't work. Error shows here:

return svd.matrixV() *  (svd.singularValues().array().abs() >
		tolerance).select(svd.singularValues().array().inverse(),
		0).matrix().asDiagonal() * svd.matrixU().adjoint();

error C2664: 'Eigen::Block<XprType,BlockRows,BlockCols,InnerPanel>::Block(const Eigen::Block<XprType,BlockRows,BlockCols,InnerPanel> &)' : cannot convert parameter 1 from 'const Eigen::GeneralProduct<Lhs,Rhs,ProductType>' to 'const Eigen::Block<XprType,BlockRows,BlockCols,InnerPanel> &'

I don't know why it occurs. Could you help me?
Comment 16 davidhigh 2014-04-14 16:55:16 UTC
In the way it is implemented now, it will probably only work correctly for matrices (but not for blocks or so). In order to get the desired behaviour, I guess one would have to follow the lines of this post: http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html#title1
Comment 17 Christoph Hertzberg 2014-04-14 17:06:56 UTC
Basically, you need to use Eigen::JacobiSVD<typename _Matrix_Type_::PlainObject>.
The return type is slightly more complicated for (partially) fixed-sized _Matrix_Type_.

But your discussion is actually not providing to fixing this bug, which would be by adding a member function to JacobiSVD (optionally also to MatrixBase).
To make this integration clean, we should wait for Bug 99 to finish (there is actually an ongoing discussion about inverse() and pseudoinverse() there).
Comment 18 Christoph Hertzberg 2014-06-20 18:41:03 UTC
The current implementation of JacobiSVD::solve appears to result in the pseudo-inverse when called with the identity matrix. We could simply define a JacobiSVD::inverse() method which returns exactly that. Any objections?
Comment 19 Jitse Niesen 2014-06-21 22:25:19 UTC
Seems reasonable.
Comment 20 Chen-Pang He 2015-02-01 08:17:37 UTC
ColPivHouseholderQR::solve returns sided inverse when called with the identity matrix.  This is the pseudo-inverse if the matrix has full rank (whether square or not).

Note You need to log in before you can comment on or make changes to this bug.