| Re: [eigen] RE: How about pseudo-inverse? |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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)
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
>>
>>
>
>
>