Hi everyone,

I just started using EIGEN and I have this code below. The purpose is to read (80,000 x 80,000) a matrix from a txt file, generate a inverse from a pseudoinverse method and write the result in a output file. I found 2 methods and I'm trying to use it ... Any comments or suggestions will be appreciate.

Thanks in advance!

Raphael

Code: Select all

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

using namespace std;
using namespace Eigen;
using Eigen::MatrixXd;
using Eigen::MatrixXf;
using Eigen::Matrix3d;

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);
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());
int main()
{
   ifstream MA;
   MA.open("A.out");

   int dimension = 8000;

   // TESTE EIGEN

   // ############################################################################################

   //template<typename Scalar>
   MatrixXf matriz;
   matriz = MatrixXf::Zero(dimension,dimension);
   //template<typename Scalar>
   MatrixXf matrizInv;
   matrizInv = 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 >> matriz(i,j);
      }
      MA.close();
   }
   else
      cout << "Unable to open file IN_AOUT";

pseudoInverse(matriz,matrizInv,3);


    return 0;
}

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::JacobiSVD< 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();
}

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::JacobiSVD< 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;
}
//-----------------------------------------------------------------------------



Hi,

for efficiency reason, I'd recommend to use a QR decomposition though the algorithm might be a little more complex.

Hi GGAEL,

Thanks for your reply. I will check this QR method.

Have a good day,

Rafael

Hi GGAEL,

I'm trying to call the function pseudoinverse() in this case, testing a small matrix, but I'm getting confused:

Code: Select all

double epslon = 2.2204e-16;
   pseudoInverse(matriz,matrizInv,epslon);

What is the correct way to call this function?

I'm also researching on the QR as you suggested.

Thanks,

Ricardo