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

[Leopold Palomo-Avellaneda <lepalom at wol dot es>]

Hi Frank,

I'm sorry for problems we have provoked you. It was my fault. The problem is 
of one of the students that we have in our lab. I simple resend his message 
because I have been in the gsl-discuss list from long time ago. I didn't know 
the existence of the gsl-help list.

Anyway, here is a new version that compiles and the same example in matlab 
code. The cpp file could be compiled in a standard linux box by:

g++ codi_gsl.cpp -lgsl -lcblas -o codi_gsl

or with gcc

gcc codi_gsl.cpp -lgsl -lcblas -lstdc++ -o codi_gsl

the file produces the result of the matrix:

0.0155688  8.90712e+17 -4.45356e+17
0.000937403 -8.50932e+17 4.25466e+17
-0.00800181 1.5173e+18 -7.58652e+17
0.0028235 -4.49257e+17 2.24628e+17
-0.000897874 -4.2349e+17 2.11745e+17
-0.000242821 1.73714e+17 -8.68568e+16
-0.000316232 -6.16453e+15 3.08227e+15
0 0 0


a similar code in matlab produces:

    0.0156    0.0012    0.0024
    0.0009    0.0070    0.0140
   -0.0080    0.0119    0.0239
    0.0028    0.0139    0.0277
   -0.0009    0.0180    0.0360
   -0.0002    0.0036    0.0072
   -0.0003    0.0000    0.0001
         0         0         0


the studend mail is attached below. 

Also, I would like to note that I have had problems with my email 
(lepalom@wol.es) so I didn't notice your mail since yesterday.


Best regards,

and thanks.

Leo


-----------------------------------------------------------------------------
Hi, good morning

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.


I need to calculate the pseudoinverse of a mxn matrix with n>m (more columns 
than rows), and with that matrix being rank-deficient.

I have created the code to compute the pseudoinverse of A. As the gsl library 
isn't able to make the SVD decomposition of n<m matrices, I used the 
transpose of A to compute the pseudoinverse of A as it follows: 

    A^T = UDV^T ==> A = VDU^T
    pinv(A) = US^(-1)V^T


The thing is that I have computed the pseudoinverse both using gsl and matlab, 
and with matlab I get more accurate results than with gsl. And my question is 
if that difference in the results is because I am coding something wrong or 
itis something else.

The result i get from matlab for the pseudoinverse of the same matrix as the 
one in the code is:

    0.0156    0.0012    0.0024
    0.0009    0.0070    0.0140
   -0.0080    0.0119    0.0239
    0.0028    0.0139    0.0277
   -0.0009    0.0180    0.0360
   -0.0002    0.0036    0.0072
   -0.0003    0.0000    0.0001
         0         0         0

I'd really apreciate your help. Thanks.

clear all;
close all;


mat =  [ 50,4.5, -23,  12,  1,  0, -1,  0;
    	 1,  2,   3,   4,  5, 1, 0, 0;
         2,  4,   6,   8, 10, 2, 0, 0];
     
% To compute the pseudoinverse in matlab I don't need to tranpose mat,
% because matlab can do the SVD decomposition of matrices with more columns
% than rows

[U,S,V] = svd(mat);


% pseudoinv = V * inv(S) * transpose(U);
n_col = length(S(:,1));
n_fil = length(S(1,:));
max = n_fil;
if (n_col < n_fil)
    max = n_col;
end

% 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
invS = transpose(S);

pseudoinv = V * invS * transpose(U)

#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
	
	gsl_matrix * SI = gsl_matrix_alloc (n_fil, n_fil);
	gsl_linalg_LU_invert (Sp, p, SI);				// Computing the inverse through LU decomposition
	gsl_permutation_free(p);
	//end Inverting S//
		
		
		
	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;
}