Pseudoinverse with svd

pers3 · **Joined:** Fri Aug 05, 2011 7:54 pm **Posts:** 1

Hello,

i implemented a method to process the pseudoinverse of a real matrix.

This is the code:

Code:

#include <iostream>

#include <stdlib.h>   //abs

#include "ap.h"
#include "ap.cpp"

#include "linalg.h"
#include "linalg.cpp"


#include "alglibinternal.h"
#include "alglibmisc.h"
#include "alglibinternal.cpp"
#include "alglibmisc.cpp"

using namespace std;
void svd(alglib::real_2d_array& a, alglib::real_2d_array& pseudoi, int m, int n)
{
      alglib::real_1d_array w;
      alglib::real_2d_array u, vt;

      alglib::rmatrixsvd(a, m, n, 2, 2, 2, w, u, vt);

      //create E_invers
      alglib::real_2d_array ei;
      ei.setlength(n, m);
      for(int i = 0; i < n; i++)
         for(int j = 0; j < m; j++)
            ei(i, j) = 0;
      for(int i = 0; i < w.length(); i++)
         if(abs(w(i)) > 10. * std::numeric_limits<double>::epsilon() )
            ei(i, i) = 1. / w(i);

      //pseudoinverse
      alglib::real_2d_array ei_ut;
      ei_ut.setlength(n, m);

      rmatrixgemm(n, m, m, 1, ei,0,0,0, u,0,0,1, 0, ei_ut,0,0);
      rmatrixgemm(n, m, n, 1, vt,0,0,1, ei_ut,0,0,0, 0, pseudoi,0,0);
}

Can somebody double-check this routine?

I tried it with

Code:

void test()
{
   int m = 5;
   int n = 5;
   alglib::real_2d_array a;
   double _r2[]   ={
               10,1,0,0,0,
               0,10,1,0,0,
               0,0,10,1,0,
               0,0,0,10,1,
               0,0,0,0,1};
   a.setcontent(m, n,_r2);

   alglib::real_2d_array pseudoi;
   pseudoi.setlength(n, m);
   svd(a, pseudoi, m, n);

   //test
   alglib::real_2d_array eye;
   eye.setlength(n, n);
   rmatrixgemm(n, n, m, 1, pseudoi,0,0,0, a,0,0,0, 0, eye,0,0);

   cout << " eye: " << endl << eye.tostring(1) << endl;
}

But the result is not convincing. Other examples are very good! It is hard to test, because some examples have a bad condition.

Thankyou for your help!!!

Erumai · **Joined:** Wed Jan 19, 2011 4:18 pm **Posts:** 11

Hi,

Sometime back I had written a small helper routine for the pseudo inverse of real matrices..and here it is..and I think it works pretty well.. I have been using it for some specific kind of image reconstructions and have even compared it with the results from matlab..they seem to agree well..

void rmatrix_pinv(alglib::real_2d_array &A,const alglib::ae_int_t m,const alglib::ae_int_t n,alglib::real_2d_array &B);
//-------------------------------------------------------------------------------
// Routine to add find the pseudoinverse of a normal matrix with linearly independent columns
// The inverse is obtained by using the LU decomposition results of ALGLIB's LU decomposition routine,
// and using the same in another routine of ALGLIB's densesolver package to solve an equation of the
// form MX=N for X, by setting N as identity matrix
// Pseudoinverse is obtained from the general formulation using inverse
// (A' . A)^-1 * A' where A' is the hermitian transpose
//
// Input Variables: A - 2D real array
// m - Number of rows in matrix A
// n - Number of columns in matrix A
// B - 2D real array to store the result
//
// On Output: B - B=pinv(A)
//
// Note: The matrix has to have linearly independent columns
//-------------------------------------------------------------------------------

void rmatrix_pinv(alglib::real_2d_array &A,const alglib::ae_int_t m,const alglib::ae_int_t n,alglib::real_2d_array &B)
{
//get the transpose and multiply it by the original matrix
alglib::real_2d_array A_tr,tmp1,tmp2;
A_tr.setlength(n,m);

tmp1.setlength(n,n);
tmp2.setlength(n,n);
alglib::rmatrixtranspose(m,n,A,0,0,A_tr,0,0); // just a transpose
int i,j;

alglib::rmatrixgemm(n,n,m,1.0,A_tr,0,0,0,A,0,0,0,0,tmp1,0,0);

alglib::ae_int_t info;
alglib::matinvreport rep;

//tmp1 has to be inverted - same as the one before
alglib::integer_1d_array pivots;
pivots.setlength(n);
rmatrixlu(tmp1,n,n,pivots);
//alglib::ae_int_t info;
alglib::densesolverreport rep_1;
alglib::real_2d_array I;
I.setlength(n,n);
eye_real(n,I);
alglib::rmatrixlusolvem(tmp1,pivots,n,I,n,info,rep_1,tmp2); //tmp2 is the inverted matrix*/

//multiply the hermitian transpose with tmp2 to get the pseudoinverse
alglib::rmatrixgemm(n,m,n,1.0,tmp2,0,0,0,A_tr,0,0,0,0,B,0,0);

}

cheers
Erumai

Alarmod · **Joined:** Wed May 23, 2012 3:06 pm **Posts:** 1

What is the function eye_real(n,I) in last topic?

forum.alglib.net

Forum rules

Pseudoinverse with svd

Who is online