Re: Errors calculating the pseudoinverse of a mxn matrix (Part II)

[Frank Reininghaus <frank78ac@googlemail.com>]

[-- Attachment #1: Type: text/plain, Size: 1358 bytes --]

Hi,

On Tuesday 08 April 2008 10:46:37 Leopold Palomo-Avellaneda wrote:
> I have been lately working with the gsl library, using the SVD
> decomposition, but I am having real trouble, so i thought that i was
> probably doing something wrong and you could help me.

yes, you are doing something wrong: You calculate the pseudoinverse in a 
completely different way in the C++ program and in the Matlab script and 
expect to get the same results. In the C++ program, you do it using 
gsl_linalg_LU_decomp () and gsl_linalg_LU_invert () which goes badly wrong 
because the third diagonal entry of the matrix S is very close to zero. In 
the Matlab script, you apparently thought of this problem and did it this 
way:
 
% As S is diagonal, que can compute its inverse by dividing inverting its
% diagonal values
for i = 1:max
   if (S(i,i) > 0.0000000001)
       S(i,i) = 1/S(i,i);
   else
       S(i,i) = 0;
   end
 end

The "if (S(i,i) > 0.0000000001)" check avoids the huge matrix entries which 
caused the problems in your C++ program. I've attached a fixed version of 
your program which reproduces the results of the Matlab script (at least up 
to the limited precision of the Matlab results).
 
Frank

P.S.: Apparently, my first reply was rejected by the list server because the 
Google Mail web interface for some reason decided to send it as text/html.

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: codi_gsl_2.cpp --]
[-- Type: text/x-c++src; charset="utf-8"; name="codi_gsl_2.cpp", Size: 3464 bytes --]

#include <cstdlib>
#include <cstdio>
#include <string>
#include <iostream>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>


int main() {

	int n_fil = 3;
	int n_col = 8;
	double mat[3][8] =  {{ 50,4.5, -23,  12,  1,  0, -1,  0},		// Example Rank-deficient matrix
		             	 {  1,  2,   3,   4,  5, 1, 0, 0},
	    	             {  2,  4,   6,   8, 10, 2, 0, 0}};
	unsigned i = 0;
	unsigned j = 0;
	gsl_matrix * gA = gsl_matrix_alloc (n_fil, n_col);
	for (i = 0; i < n_fil; i++)
		for (j = 0; j < n_col; j++)
	   		gsl_matrix_set (gA, i, j, mat[i][j]);


	gsl_matrix * gA_t = gsl_matrix_alloc (n_col, n_fil);
	gsl_matrix_transpose_memcpy (gA_t, gA);					// Computing the transpose of gA
		
	gsl_matrix * U = gsl_matrix_alloc (n_col, n_fil);
	gsl_matrix * V= gsl_matrix_alloc (n_fil, n_fil);
	gsl_vector * S = gsl_vector_alloc (n_fil);


	// Computing the SVD of the transpose of A
	// The matrix 'gA_t' will contain 'U' after the function is called
	gsl_vector * work = gsl_vector_alloc (n_fil);
	gsl_linalg_SV_decomp (gA_t, V, S, work);
	gsl_vector_free(work);
		
	gsl_matrix_memcpy (U, gA_t);
	
				
	//Inverting S//
	//----------------------------------------------------------
	// Matrix 'S' is diagonal, so it is contained in a vector.
	// We operate to convert the vector 'S' into the matrix 'Sp'.
	//Then we invert 'Sp' to 'Spu'
	//----------------------------------------------------------
	gsl_matrix * Sp = gsl_matrix_alloc (n_fil, n_fil);
	gsl_matrix_set_zero (Sp);
	for (i = 0; i < n_fil; i++)
		gsl_matrix_set (Sp, i, i, gsl_vector_get(S, i));	// Vector 'S' to matrix 'Sp'
	
	gsl_permutation * p = gsl_permutation_alloc (n_fil);
	int signum;
	gsl_linalg_LU_decomp (Sp, p, &signum);				// Computing the LU decomposition

	// Compute the inverse like in the MATLAB script

	gsl_matrix * SI = gsl_matrix_calloc (n_fil, n_fil);

	for (i = 0; i < n_fil; i++) {
	  std::cout << "S [" << i << "] = " << gsl_vector_get (S, i) << std::endl;

	  if (gsl_vector_get (S, i) > 0.0000000001)
	    gsl_matrix_set (SI, i, i, 1.0 / gsl_vector_get (S, i));
	}		
	
	gsl_matrix * VT = gsl_matrix_alloc (n_fil, n_fil);
	gsl_matrix_transpose_memcpy (VT, V);					// Tranpose of V
		
		
	//THE PSEUDOINVERSE//
	//----------------------------------------------------------
	//Computation of the pseudoinverse of trans(A) as pinv(A) = U·inv(S).trans(V)   with trans(A) = U.S.trans(V)
	//----------------------------------------------------------
	gsl_matrix * SIpVT = gsl_matrix_alloc (n_fil, n_fil);
	gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,				// Calculating  inv(S).trans(V)
                	1.0, SI, VT,
                	0.0, SIpVT);

			
	gsl_matrix * pinv = gsl_matrix_alloc (n_col, n_fil);	// Calculating  U·inv(S).trans(V)
	gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                	1.0, U, SIpVT,
                	0.0, pinv);

	//end THE PSEUDOINVERSE//
	
   	std::cout << "pinv:" << std::endl;
   	for (i = 0; i < n_col; i++)  
    		for (j = 0; j < n_fil; j++)
    			printf ("m(%d,%d) = %g\n", i, j, 
    				gsl_matrix_get (pinv, i, j));
   	std::cout << "\n" << std::endl;
	


	gsl_matrix_free(VT);
	gsl_matrix_free(SI);
	gsl_matrix_free(SIpVT);
	gsl_matrix_free(gA_t);
	gsl_matrix_free(U);
	gsl_matrix_free(gA);
	gsl_matrix_free(V);
	gsl_vector_free(S);




	return 0;
}