From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 20015 invoked by alias); 31 Mar 2008 18:36:19 -0000 Received: (qmail 20006 invoked by uid 22791); 31 Mar 2008 18:36:18 -0000 X-Spam-Check-By: sourceware.org Received: from aybabtu.com (HELO aybabtu.com) (69.60.117.155) by sourceware.org (qpsmtpd/0.31) with ESMTP; Mon, 31 Mar 2008 18:35:58 +0000 Received: from 62.57.84.101.dyn.user.ono.com ([62.57.84.101] helo=cami.somiatruites.net) by aybabtu.com with esmtps (TLS-1.0:RSA_AES_256_CBC_SHA1:32) (Exim 4.63) (envelope-from ) id 1JgOrG-0007i3-If for gsl-discuss@sourceware.org; Mon, 31 Mar 2008 20:35:51 +0200 Received: from indiana.somiatruites.net ([192.168.1.8] ident=leo) by cami.somiatruites.net with esmtp (Exim 4.63) (envelope-from ) id 1JgOqd-00020v-Tn for gsl-discuss@sourceware.org; Mon, 31 Mar 2008 20:35:11 +0200 From: Leopold Palomo Avellaneda Reply-To: lepalom@wol.es To: GSL Discuss Mailing List Subject: Errors calculating the pseudoinverse of a mxn matrix Date: Mon, 31 Mar 2008 18:36:00 -0000 User-Agent: KMail/1.9.7 MIME-Version: 1.0 Content-Type: multipart/signed; boundary="nextPart1276729.D2AuRmdJ2g"; protocol="application/pgp-signature"; micalg=pgp-sha1 Content-Transfer-Encoding: 7bit Message-Id: <200803312035.10232.lepalom@wol.es> 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-q1/txt/msg00047.txt.bz2 --nextPart1276729.D2AuRmdJ2g Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: quoted-printable Content-Disposition: inline Content-length: 4012 Hi, good morning I have been lately working with the gsl library, using the SVD decompositio= n,=20 but I am having real trouble, so i thought that i was probably doing=20 something wrong and you could help me. I need to calculate the pseudoinverse of a mxn matrix with n>m (more column= s=20 than rows), and with that matrix being rank-deficient. I have created the code, that you can read next, to compute the pseudoinver= se=20 of A. As the gsl library isn't able to make the SVD decomposition of n A =3D VDU^T pinv(A) =3D US^(-1)V^T The pseudoinverse is correct when the A is full-rank, but not when is=20 rank-deficient. And I don't know what i am doing wrong. The code I've written is: int n_fil =3D 3; int n_col =3D 8; double mat[3][8] =3D {{ 50,4.5, -23, 12, 1, 0, -1, 0}, // Examp= le=20 Rank-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= =20 transpose of gA =20=20=20=20=20=20=20=20 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); =20=20=20=20=20=20=20=20 U =3DgA_t; =20=20=20=20 =20=20=20=20=20=20=20=20=20=20=20=20=20=20=20=20 //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=20 matrix 'Sp' =20=20=20=20 gsl_permutation * p =3D gsl_permutation_alloc (n_fil); int signum; gsl_linalg_LU_decomp (Sp, p, &signum); // Computing the LU=20 decomposition gsl_matrix * SI =3D gsl_matrix_alloc (n_fil, n_fil); gsl_linalg_LU_invert (Sp, p, SI); // Computing the inverse=20 through LU decomposition gsl_permutation_free(p); //end Inverting S// =20=20=20=20=20=20=20=20 =20=20=20=20=20=20=20=20 =20=20=20=20=20=20=20=20 gsl_matrix * VT =3D gsl_matrix_alloc (n_fil, n_fil); gsl_matrix_transpose_memcpy (VT, V); // Tranpose of V =20=20=20=20=20=20=20=20 =20=20=20=20=20=20=20=20 //THE PSEUDOINVERSE// //---------------------------------------------------------- //Computation of the pseudoinverse of trans(A) as pinv(A) =3D U=C2=B7inv(S)= .trans(V)=20=20=20 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= =20=20 inv(S).trans(V) 1.0, SI, VT, 0.0, SIpVT); =20=20=20=20=20=20=20=20 gsl_matrix * pinv =3D gsl_matrix_alloc (n_col, n_fil); // Calculating=20= =20 U=C2=B7inv(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// --=20 -- Linux User 152692 PGP: 0xF944807E Catalonia --nextPart1276729.D2AuRmdJ2g Content-Type: application/pgp-signature; name=signature.asc Content-Description: This is a digitally signed message part. Content-length: 189 -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.6 (GNU/Linux) iD8DBQBH8S7eaC7AnPlEgH4RAl+LAJ9vkL0hN9VzF7UcO9q8F2vSpOSddQCgouqj eXW3zwbaT9O91HTxjgbrZUE= =v1VK -----END PGP SIGNATURE----- --nextPart1276729.D2AuRmdJ2g--