From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 6510 invoked by alias); 10 Apr 2008 19:17:00 -0000 Received: (qmail 6499 invoked by uid 22791); 10 Apr 2008 19:16:58 -0000 X-Spam-Check-By: sourceware.org Received: from mu-out-0910.google.com (HELO mu-out-0910.google.com) (209.85.134.184) by sourceware.org (qpsmtpd/0.31) with ESMTP; Thu, 10 Apr 2008 19:16:21 +0000 Received: by mu-out-0910.google.com with SMTP id g7so197177muf.0 for ; Thu, 10 Apr 2008 12:16:17 -0700 (PDT) Received: by 10.82.154.5 with SMTP id b5mr2881077bue.10.1207854977217; Thu, 10 Apr 2008 12:16:17 -0700 (PDT) Received: from noname ( [87.189.11.199]) by mx.google.com with ESMTPS id j2sm6351324mue.3.2008.04.10.12.16.13 (version=TLSv1/SSLv3 cipher=OTHER); Thu, 10 Apr 2008 12:16:16 -0700 (PDT) From: Frank Reininghaus To: Leopold Palomo-Avellaneda Subject: Re: Errors calculating the pseudoinverse of a mxn matrix (Part II) Date: Thu, 10 Apr 2008 19:17:00 -0000 User-Agent: KMail/1.9.6 (enterprise 0.20070907.709405) Cc: gsl-discuss@sourceware.org, leo@alaxarxa.net, Josep Arnau References: <200804081046.37469.lepalom@wol.es> In-Reply-To: <200804081046.37469.lepalom@wol.es> MIME-Version: 1.0 X-Length: 6397 X-UID: 34 Content-Type: Multipart/Mixed; boundary="Boundary-00=_8dm/HJOnQhGBMiC" Message-Id: <200804102116.12005.frank78ac@googlemail.com> Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q2/txt/msg00013.txt.bz2 --Boundary-00=_8dm/HJOnQhGBMiC Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: 7bit Content-Disposition: inline Content-length: 1358 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. --Boundary-00=_8dm/HJOnQhGBMiC Content-Type: text/x-c++src; charset="utf-8"; name="codi_gsl_2.cpp" Content-Transfer-Encoding: quoted-printable Content-Disposition: attachment; filename="codi_gsl_2.cpp" Content-length: 3473 #include #include #include #include #include #include #include #include int main() { int n_fil =3D 3; int n_col =3D 8; double mat[3][8] =3D {{ 50,4.5, -23, 12, 1, 0, -1, 0}, // Example Ra= nk-deficient matrix { 1, 2, 3, 4, 5, 1, 0, 0}, { 2, 4, 6, 8, 10, 2, 0, 0}}; unsigned i =3D 0; unsigned j =3D 0; gsl_matrix * gA =3D gsl_matrix_alloc (n_fil, n_col); for (i =3D 0; i < n_fil; i++) for (j =3D 0; j < n_col; j++) gsl_matrix_set (gA, i, j, mat[i][j]); gsl_matrix * gA_t =3D gsl_matrix_alloc (n_col, n_fil); gsl_matrix_transpose_memcpy (gA_t, gA); // Computing the transpose of = gA =09=09 gsl_matrix * U =3D gsl_matrix_alloc (n_col, n_fil); gsl_matrix * V=3D gsl_matrix_alloc (n_fil, n_fil); gsl_vector * S =3D 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 =3D gsl_vector_alloc (n_fil); gsl_linalg_SV_decomp (gA_t, V, S, work); gsl_vector_free(work); =09=09 gsl_matrix_memcpy (U, gA_t); =09 =09=09=09=09 //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 =3D gsl_matrix_alloc (n_fil, n_fil); gsl_matrix_set_zero (Sp); for (i =3D 0; i < n_fil; i++) gsl_matrix_set (Sp, i, i, gsl_vector_get(S, i)); // Vector 'S' to matrix = 'Sp' =09 gsl_permutation * p =3D 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 =3D gsl_matrix_calloc (n_fil, n_fil); for (i =3D 0; i < n_fil; i++) { std::cout << "S [" << i << "] =3D " << gsl_vector_get (S, i) << std::end= l; if (gsl_vector_get (S, i) > 0.0000000001) gsl_matrix_set (SI, i, i, 1.0 / gsl_vector_get (S, i)); }=09=09 =09 gsl_matrix * VT =3D gsl_matrix_alloc (n_fil, n_fil); gsl_matrix_transpose_memcpy (VT, V); // Tranpose of V =09=09 =09=09 //THE PSEUDOINVERSE// //---------------------------------------------------------- //Computation of the pseudoinverse of trans(A) as pinv(A) =3D U=B7inv(S).t= rans(V) with trans(A) =3D U.S.trans(V) //---------------------------------------------------------- gsl_matrix * SIpVT =3D gsl_matrix_alloc (n_fil, n_fil); gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, // Calculating inv(S).tran= s(V) 1.0, SI, VT, 0.0, SIpVT); =09=09=09 gsl_matrix * pinv =3D gsl_matrix_alloc (n_col, n_fil); // Calculating U= =B7inv(S).trans(V) gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, U, SIpVT, 0.0, pinv); //end THE PSEUDOINVERSE// =09 std::cout << "pinv:" << std::endl; for (i =3D 0; i < n_col; i++)=20=20 for (j =3D 0; j < n_fil; j++) printf ("m(%d,%d) =3D %g\n", i, j,=20 gsl_matrix_get (pinv, i, j)); std::cout << "\n" << std::endl; =09 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; } --Boundary-00=_8dm/HJOnQhGBMiC--