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
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)
Parts of this can ideally be used to resolve bug 241 (using LDLT to solve linear equations with singular matrices).
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(); }
(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);
My own pseudoinverse looks like this and yes, this code will work only with ThinU and ThinV.
*** Bug 412 has been marked as a duplicate of this bug. ***
(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.
(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.
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.
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
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-------------------------------------------------
(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.
Thanks Christoph for the hint! It works now.
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
(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?
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
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).
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?
Seems reasonable.
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).
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!
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
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
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; }
BTW, you can also directly call: MatrixXd A = ...; MatrixXd pinv = A.completeOrthogonalDecomposition().pseudoInverse(); but so far this path also requires Dynamic-sized matrices.
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
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
-- 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.