RE: [eigen] How about pseudo-inverse?

[Yohann Solaro <ysolaro@xxxxxxxxx>]

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