Errors calculating the pseudoinverse of a mxn matrix

Errors calculating the pseudoinverse of a mxn matrix

[Leopold Palomo Avellaneda]

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

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, that you can read next, 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 pseudoinverse is correct when the A is full-rank, but not when is 
rank-deficient. And I don't know what i am doing wrong.


The code I've written is:



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);
        
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);


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);

//end THE PSEUDOINVERSE//
-- 
--
Linux User 152692
PGP: 0xF944807E
Catalonia

[-- Attachment #2: This is a digitally signed message part. --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

Re: Errors calculating the pseudoinverse of a mxn matrix

[Frank Reininghaus]

Hi Leopold,

when asking questions about unexpected behaviour of a program, you should 
always send the *complete* and *compilable* source code. Looking up the 
header files that are required by your program and removing all the comments 
with incorrect line breaks that cause compiler errors isn't much fun for 
people who are willing to help you. Moreover, it would be helpful if the 
program produced some output that clearly shows the unexpected result.

BTW: Questions like this should be sent to the 'gsl-help' mailing list:

http://lists.gnu.org/mailman/listinfo/help-gsl

Regards,
Frank

On Monday 31 March 2008 20:35:06 Leopold Palomo Avellaneda wrote:
> 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.