# Errors calculating the pseudoinverse of a mxn matrix

Leopold Palomo Avellaneda lepalom@wol.es
Mon Mar 31 18:36:00 GMT 2008

```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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part.
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20080331/c448cf10/attachment.sig>
```