Re: [eigen] RE: How about pseudo-inverse?

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/1/22 Andre Krause <post@xxxxxxxxxxxxxxxx>:
> what is good about forcing "Eigen::RowMajor", when the eigen default is column
> major ? btw, is there an easy function that allows me to convert a matrix from
> column to row major ? (i just use not that often, and didnt found that info in
> the docs at first glance)

You can simply do:

ColMajorMatrixType m;
......
RowMajorMatrixType m2 = m;

Benoit


>
> hamelin.philippe@xxxxxxx wrote:
>> Ok now we have two versions: the one from Yohann and the one from Ricard.. I would like to suggest a new function with I think the best of both. Firstly, here are my comments:
>>
>> Ricard:
>>
>> * Pros:
>>  - Your function is more compact (elegant)
>> * Cons:
>>  - You don't take the absolute value of the eigen values for the tolerance check
>>
>> Yohann:
>>
>> * Pros:
>>  - You take the absolute values for the tolerance check
>>  - You force Eigen::RowMajor which I think is a good practice
>> * Cons:
>>  - Your function is less compact
>>
>> So here is the combination of both functions:
>>
>>     template<typename Scalar>
>>     static bool pseudoInverse(const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> &a,
>>         Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> &result,
>>         double epsilon = std::numeric_limits<Scalar>::epsilon())
>>     {
>>         if(a.rows()<a.cols())
>>             return false;
>>
>>         Eigen::SVD< Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > svd = a.svd();
>>
>>         Scalar tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().cwise().abs().maxCoeff();
>>
>>         result = svd.matrixV() * (svd.singularValues().cwise()..abs() > tolerance).select(svd.singularValues().
>>             cwise().inverse(), 0).asDiagonal() * svd.matrixU().adjoint();
>>     }
>>
>> I used std::numeric_limits<>::epsilon() to have a default value for epsilon. Compiled but not yet tested. Your comments are very welcome!
>>
>> Philippe
>>
>> -----Message d'origine-----
>> De : Listengine [mailto:listengine@xxxxxxxxxxxxxxxxx] De la part de Yohann Solaro
>> Envoyé : 22 janvier 2010 04:49
>> À : eigen@xxxxxxxxxxxxxxxxxxx
>> Objet : RE: [eigen] How about pseudo-inverse?
>>
>> Hi,
>>
>> We are using Eigen at Movea, since few months, for Motion analysis .
>> We needed the pseudo-inverse computation and decided to implement it using the SVD module.
>> The code is not really compact and I think we can do more optimizations.
>> There is a bit of optimization on "mAdjointU" that consists in shortening its size because we noticed that it contains a lot of zero (sparse ?).
>> Our main source for the theory was : http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#The_general_case_and_the_SVD_method
>>
>>
>> //-----------------------------------------------------------------------------
>> //! \brief pseudo_inverse
>> //! \brief Computes the pseudo-Inverse of a rectangular Matrix
>> //! \return true if the the SVD succeed
>> //! \param[in] The Original Matrix
>> //! \param[in] The Pseudo-Inverted Matrix
>> //! \tparam The type contained by the matrix (double, int...)
>> template<typename Scalar>
>> bool pinv(const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> &a, Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> &a_pinv)
>> {
>>     // see : http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#The_general_case_and_the_SVD_method
>>
>>     if ( a.rows()<a.cols() )
>>         return false;
>>
>>     // SVD
>>     Eigen::SVD< Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> > svdA(a);
>>
>>     Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> vSingular = svdA.singularValues();
>>
>>     // Build a diagonal matrix with the Inverted Singular values
>>     // The pseudo inverted singular matrix is easy to compute :
>>     // is formed by replacing every nonzero entry by its reciprocal (inversing).
>>     Eigen::Matrix<Scalar, Eigen::Dynamic, 1, Eigen::RowMajor> vPseudoInvertedSingular(svdA.matrixV().cols(),1);
>>
>>     for (int iRow =0; iRow<vSingular.rows(); iRow++)
>>     {
>>         if ( fabs(vSingular(iRow))<=1e-10 ) // Todo : Put epsilon in parameter
>>         {
>>             vPseudoInvertedSingular(iRow,0)=0.;
>>         }
>>         else
>>         {
>>             vPseudoInvertedSingular(iRow,0)=1./vSingular(iRow);
>>         }
>>     }
>>
>>     // A little optimization here
>>     Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mAdjointU = svdA.matrixU().adjoint().block(0,0,vSingular.rows(),svdA.matrixU().adjoint().cols());
>>
>>     // Pseudo-Inversion : V * S * U'
>>     a_pinv = (svdA.matrixV() *  vPseudoInvertedSingular.asDiagonal()) * mAdjointU  ;
>>
>>     return true;
>> }
>> //-----------------------------------------------------------------------------
>>
>>
>>
>>
>> Feel free to use it anywhere you want,
>> Feedback are welcome.
>>
>>
>>
>> Best Regards,
>>
>>
>>
>>
>> Yohann Solaro
>> Apprentice Software Engineer
>>
>> Movea SA
>> 7, Parvis Louis Néel, BP50 F-38040 Grenoble
>> Email:     ysolaro@xxxxxxxxx
>> Web:       www.movea.com
>>                 www.gyration.com
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> -----Message d'origine-----
>> De : Listengine [mailto:listengine@xxxxxxxxxxxxxxxxx] De la part de Ricard Marxer Piñón
>> Envoyé : jeudi 21 janvier 2010 16:16
>> À : eigen@xxxxxxxxxxxxxxxxxxx
>> Objet : Re: [eigen] How about pseudo-inverse?
>>
>> Hi,
>>
>> I implemented my own version like this.  But I don't know if it is
>> really correct, so please correct me if I'm wrong:
>>
>> #include <Eigen/Core>
>> #include <Eigen/Array>
>> #include <Eigen/SVD>
>>
>> typedef Eigen::Matrix< Real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR;
>>
>> // SVD method
>> void pseudoInverse(const MatrixXR& a, MatrixXR* result, Real epsilon) {
>>   Eigen::SVD<MatrixXR> svd = a.svd();
>>
>>   // As done in Matlab:
>>   // http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
>>   Real tolerance = epsilon * max(a.cols(), a.rows()) *
>> svd.singularValues().cwise().abs().maxCoeff();
>>
>>   (*result) = svd.matrixV() * (svd.singularValues().cwise() >
>> tolerance).select(svd.singularValues().cwise().inverse(),
>> 0).asDiagonal() * svd.matrixU().adjoint();
>> }
>>
>> ricard
>>
>> On Thu, Jan 21, 2010 at 3:55 PM,  <hamelin.philippe@xxxxxxx> wrote:
>>> Hello,
>>> I'm just new to the mailing list (and to Eigen2). I'm currently programming
>>> an underwater vehicle simulator using Eigen2. I actually need the
>>> pseudo-inverse of a rectangular matrix. I didn't find anything about that in
>>> the code. There's seems to be a preliminary version of the SVD which could
>>> help to do pseudo-inverse. So, does anyone has something to suggest to me
>>> about that?
>>>
>>> Thank you,
>>>
>>> Philippe Hamelin, ing. jr, M. Ing
>>> Chercheur / Researcher
>>>
>>> Tel: (450) 652-8499 x2198
>>> Fax: (450) 652-1316
>>>
>>> Institut de recherche d'Hydro-Québec
>>> Unité robotique et civil
>>> 1740, boul. Lionel-Boulet
>>> Varennes (QC) J1X 1S1, Canada
>>>
>>>
>>
>>
>>
>
>
>
>