This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 257 - pseudoinverse
Summary: pseudoinverse
Status: CONFIRMED
Alias: None
Product: Eigen
Classification: Unclassified
Component: SVD (show other bugs)
Version: 3.0
Hardware: All All
: Low Feature Request
Assignee: Nobody
URL:
Whiteboard:
Keywords:
: 412 (view as bug list)
Depends on: ExprEval
Blocks:
  Show dependency treegraph
 
Reported: 2011-04-28 15:52 UTC by Cyrille Berger
Modified: 2019-12-04 10:41 UTC (History)
16 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).
Comment 21 Joseph Dien 2019-06-21 18:32:13 UTC
Hi, I would just like to bump up this request now that Bug 99 has apparently been resolved.  The code snippet in the FAQ is not at all easy to use for a newbie such as myself.  Thanks for such an amazing resource!
Comment 22 Sylvain Rousseau 2019-11-09 17:13:44 UTC
Hi !

I tried to compile the snippet given @ comment14, but it fails.

The error is related to the result matrix type which cannot be the same as the input matrix.
Indeed, the dimensions of the pseudo inverse are switched with respect to the input matrix ones.

I use Eigen dev version.

May I ask you some help ?

Below the compiler error message :

                                                   ^
In file included from /usr/local/include/Eigen/Core:152:0,
                 from ../../CPhaseSensor/include/PhaseSensorDefs.h:28,
                 from ../include/COpdCtrlDefs.h:24,
                 from ../include/COpdCtrl.h:25,
                 from ../src/COpdCtrl.cpp:35:
/usr/local/include/Eigen/src/Core/GeneralProduct.h: In instantiation of ‘const Eigen::Product<Derived, OtherDerived> Eigen::MatrixBase<Derived>::operator*(const Eigen::MatrixBase<OtherDerived>&) const [with OtherDerived = Eigen::Transpose<const Eigen::Matrix<double, 15, 15, 0, 15, 15> >; Derived = Eigen::Product<Eigen::Matrix<double, 6, 6, 0, 6, 6>, Eigen::DiagonalWrapper<const Eigen::MatrixWrapper<const Eigen::Select<Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<double, double, (Eigen::internal::ComparisonName)1u>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Array<double, 6, 1, 0, 6, 1> >, const Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, 6, 1, 0, 6, 1> > > >, Eigen::CwiseUnaryOp<Eigen::internal::scalar_inverse_op<double>, const Eigen::ArrayWrapper<const Eigen::Matrix<double, 6, 1, 0, 6, 1> > >, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Array<double, 6, 1, 0, 6, 1> > > > >, 1>]’:
../../CPhaseSensor/include/PhaseSensorDefs.h:261:2:   required from ‘MatType pseudoInverse(const MatType&, double) [with MatType = Eigen::Matrix<double, 15, 6>]’
../src/COpdCtrl.cpp:81:31:   required from here
/usr/local/include/Eigen/src/Core/GeneralProduct.h:421:3: error: static assertion failed: INVALID_MATRIX_PRODUCT
   EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)

Cheers!
Sylvain
Comment 23 Sylvain Rousseau 2019-11-09 17:58:59 UTC
To complete my previous comment : it sounds like the pseudoInverse() function code proposed is not compliant with the statically sized matrices which explains the compiler compliance.

The underlying reason seems to be that Eigen::JacobiSVD object does not accept static matrices when called with the ComputeThinU | ComputeThinV flags.

Could you confirm ?

Cheers
Comment 24 Gael Guennebaud 2019-11-13 15:13:35 UTC
Here is an updated version improved in several ways:

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


template<typename MatType>
using PseudoInverseType = Eigen::Matrix<typename MatType::Scalar, MatType::ColsAtCompileTime, MatType::RowsAtCompileTime>;

template<typename MatType>
PseudoInverseType<MatType> pseudoInverse(const MatType &a, double epsilon = std::numeric_limits<double>::epsilon())
{
	using WorkingMatType = Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, Eigen::Dynamic, 0,
																			 MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime>;
	Eigen::BDCSVD<WorkingMatType> svd(a, Eigen::ComputeThinU | Eigen::ComputeThinV);
	svd.setThreshold(epsilon*std::max(a.cols(), a.rows()));
	Eigen::Index rank = svd.rank();
	Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, MatType::RowsAtCompileTime,
								0, Eigen::BDCSVD<WorkingMatType>::MaxDiagSizeAtCompileTime, MatType::MaxRowsAtCompileTime>
		tmp = svd.matrixU().leftCols(rank).adjoint();
	tmp = svd.singularValues().head(rank).asDiagonal().inverse() * tmp;
	return svd.matrixV().leftCols(rank) * tmp;
}


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

	Eigen::Matrix<double,3,2> 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;
}
Comment 25 Gael Guennebaud 2019-11-13 15:18:16 UTC
BTW, you can also directly call:

    MatrixXd A = ...;
    MatrixXd pinv = A.completeOrthogonalDecomposition().pseudoInverse();

but so far this path also requires Dynamic-sized matrices.
Comment 26 Sylvain Rousseau 2019-11-15 13:54:57 UTC
Hi Gael,

Thank you very much for your answer.

I prefere to avoid completeOrthogonalDecomposition() which is known for its higher complexity than SVD algorithm.

Cheers !

Sylvain
Comment 27 Sylvain Rousseau 2019-11-18 11:18:27 UTC
Hi Gael,

I'd have two questions about your code :

- 1 -
using WorkingMatType = Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, Eigen::Dynamic, 0, 
MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime>;

I've been wondering what's the purpose of the options "MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime" ?

- 2 -
Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, MatType::RowsAtCompileTime, 0, Eigen::JacobiSVD<WorkingMatType>::MaxDiagSizeAtCompileTime, MatType::MaxRowsAtCompileTime>

Why do you define this particular type ? 
What's the advantage compared to a simple dynamic matrix type like Matrix<double, Dynamic, Dynamic > ?

Thanks for your explanations.

Cheers

Sylvain
Comment 28 Nobody 2019-12-04 10:41:40 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/257.

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