* Errors calculating the pseudoinverse of a mxn matrix (Part II)
@ 2008-04-08 8:47 Leopold Palomo-Avellaneda
2008-04-10 19:17 ` Frank Reininghaus
0 siblings, 1 reply; 2+ messages in thread
From: Leopold Palomo-Avellaneda @ 2008-04-08 8:47 UTC (permalink / raw)
To: gsl-discuss, lepalom; +Cc: Frank Reininghaus, leo, Josep Arnau
[-- Attachment #1: Type: text/plain, Size: 2673 bytes --]
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.
[-- Attachment #2: matlab_code.m --]
[-- Type: text/x-objcsrc, Size: 731 bytes --]
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)
[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #3: codi_gsl.cpp --]
[-- Type: text/x-c++src; charset="utf-8"; name="codi_gsl.cpp", Size: 3329 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
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;
}
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: Errors calculating the pseudoinverse of a mxn matrix (Part II)
2008-04-08 8:47 Errors calculating the pseudoinverse of a mxn matrix (Part II) Leopold Palomo-Avellaneda
@ 2008-04-10 19:17 ` Frank Reininghaus
0 siblings, 0 replies; 2+ messages in thread
From: Frank Reininghaus @ 2008-04-10 19:17 UTC (permalink / raw)
To: Leopold Palomo-Avellaneda; +Cc: gsl-discuss, leo, Josep Arnau
[-- 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;
}
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2008-04-10 19:17 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2008-04-08 8:47 Errors calculating the pseudoinverse of a mxn matrix (Part II) Leopold Palomo-Avellaneda
2008-04-10 19:17 ` Frank Reininghaus
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).