[eigen] RE: How about pseudo-inverse?

[<hamelin.philippe@xxxxxxx>]

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
>
>



-- 
ricard
http://www.ricardmarxer.com
http://www.caligraft.com