From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 31403 invoked by alias); 8 Apr 2008 08:47:14 -0000 Received: (qmail 31375 invoked by uid 22791); 8 Apr 2008 08:47:12 -0000 X-Spam-Check-By: sourceware.org Received: from dash.upc.es (HELO dash.upc.es) (147.83.2.50) by sourceware.org (qpsmtpd/0.31) with ESMTP; Tue, 08 Apr 2008 08:46:51 +0000 Received: from hamilton.upcnetadm.upcnet.es (hamilton.upcnetadm.upcnet.es [147.83.2.240]) by dash.upc.es (8.14.1/8.13.1) with ESMTP id m388kbMP014475; Tue, 8 Apr 2008 10:46:38 +0200 Received: from jargon.upc.es ([147.83.37.191]) by hamilton.upcnetadm.upcnet.es (Lotus Domino Release 5.0.12) with ESMTP id 2008040810463814:1764 ; Tue, 8 Apr 2008 10:46:38 +0200 From: Leopold Palomo-Avellaneda To: gsl-discuss@sourceware.org, lepalom@wol.es Subject: Errors calculating the pseudoinverse of a mxn matrix (Part II) Date: Tue, 08 Apr 2008 08:47:00 -0000 User-Agent: KMail/1.9.7 Cc: Frank Reininghaus , leo@alaxarxa.net, Josep Arnau MIME-Version: 1.0 Message-Id: <200804081046.37469.lepalom@wol.es> Content-Type: Multipart/Mixed; boundary="Boundary-00=_tDz+HNGlP3bHTH2" X-Mail-Scanned: Criba 2.0 + Clamd X-Greylist: Sender IP whitelisted, not delayed by milter-greylist-3.0 (dash.upc.es [147.83.2.50]); Tue, 08 Apr 2008 10:46:38 +0200 (CEST) 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/msg00007.txt.bz2 --Boundary-00=_tDz+HNGlP3bHTH2 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset="utf-8" Content-Disposition: inline Content-length: 2673 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 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. --Boundary-00=_tDz+HNGlP3bHTH2 Content-Transfer-Encoding: 7bit Content-Type: text/x-objcsrc; charset="utf-8"; name="matlab_code.m" Content-Disposition: attachment; filename="matlab_code.m" Content-length: 731 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) --Boundary-00=_tDz+HNGlP3bHTH2 Content-Type: text/x-c++src; charset="utf-8"; name="codi_gsl.cpp" Content-Disposition: attachment; filename="codi_gsl.cpp" Content-Transfer-Encoding: base64 Content-length: 4514 I2luY2x1ZGUgPGNzdGRsaWI+DQojaW5jbHVkZSA8Y3N0ZGlvPg0KI2luY2x1 ZGUgPHN0cmluZz4NCiNpbmNsdWRlIDxpb3N0cmVhbT4NCg0KI2luY2x1ZGUg PGdzbC9nc2xfbWF0cml4Lmg+DQojaW5jbHVkZSA8Z3NsL2dzbF9wZXJtdXRh dGlvbi5oPg0KI2luY2x1ZGUgPGdzbC9nc2xfbGluYWxnLmg+DQojaW5jbHVk ZSA8Z3NsL2dzbF9ibGFzLmg+DQoNCg0KaW50IG1haW4oKSB7DQoNCglpbnQg bl9maWwgPSAzOw0KCWludCBuX2NvbCA9IDg7DQoJZG91YmxlIG1hdFszXVs4 XSA9ICB7eyA1MCw0LjUsIC0yMywgIDEyLCAgMSwgIDAsIC0xLCAgMH0sCQkv LyBFeGFtcGxlIFJhbmstZGVmaWNpZW50IG1hdHJpeA0KCQkgICAgICAgICAg ICAgCSB7ICAxLCAgMiwgICAzLCAgIDQsICA1LCAxLCAwLCAwfSwNCgkgICAg CSAgICAgICAgICAgICB7ICAyLCAgNCwgICA2LCAgIDgsIDEwLCAyLCAwLCAw fX07DQoJdW5zaWduZWQgaSA9IDA7DQoJdW5zaWduZWQgaiA9IDA7DQoJZ3Ns X21hdHJpeCAqIGdBID0gZ3NsX21hdHJpeF9hbGxvYyAobl9maWwsIG5fY29s KTsNCglmb3IgKGkgPSAwOyBpIDwgbl9maWw7IGkrKykNCgkJZm9yIChqID0g MDsgaiA8IG5fY29sOyBqKyspDQoJICAgCQlnc2xfbWF0cml4X3NldCAoZ0Es IGksIGosIG1hdFtpXVtqXSk7DQoNCg0KCWdzbF9tYXRyaXggKiBnQV90ID0g Z3NsX21hdHJpeF9hbGxvYyAobl9jb2wsIG5fZmlsKTsNCglnc2xfbWF0cml4 X3RyYW5zcG9zZV9tZW1jcHkgKGdBX3QsIGdBKTsJCQkJCS8vIENvbXB1dGlu ZyB0aGUgdHJhbnNwb3NlIG9mIGdBDQoJCQ0KCWdzbF9tYXRyaXggKiBVID0g Z3NsX21hdHJpeF9hbGxvYyAobl9jb2wsIG5fZmlsKTsNCglnc2xfbWF0cml4 ICogVj0gZ3NsX21hdHJpeF9hbGxvYyAobl9maWwsIG5fZmlsKTsNCglnc2xf dmVjdG9yICogUyA9IGdzbF92ZWN0b3JfYWxsb2MgKG5fZmlsKTsNCg0KDQoJ Ly8gQ29tcHV0aW5nIHRoZSBTVkQgb2YgdGhlIHRyYW5zcG9zZSBvZiBBDQoJ Ly8gVGhlIG1hdHJpeCAnZ0FfdCcgd2lsbCBjb250YWluICdVJyBhZnRlciB0 aGUgZnVuY3Rpb24gaXMgY2FsbGVkDQoJZ3NsX3ZlY3RvciAqIHdvcmsgPSBn c2xfdmVjdG9yX2FsbG9jIChuX2ZpbCk7DQoJZ3NsX2xpbmFsZ19TVl9kZWNv bXAgKGdBX3QsIFYsIFMsIHdvcmspOw0KCWdzbF92ZWN0b3JfZnJlZSh3b3Jr KTsNCgkJDQoJZ3NsX21hdHJpeF9tZW1jcHkgKFUsIGdBX3QpOw0KCQ0KCQkJ CQ0KCS8vSW52ZXJ0aW5nIFMvLw0KCS8vLS0tLS0tLS0tLS0tLS0tLS0tLS0t LS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQ0KCS8vIE1h dHJpeCAnUycgaXMgZGlhZ29uYWwsIHNvIGl0IGlzIGNvbnRhaW5lZCBpbiBh IHZlY3Rvci4NCgkvLyBXZSBvcGVyYXRlIHRvIGNvbnZlcnQgdGhlIHZlY3Rv ciAnUycgaW50byB0aGUgbWF0cml4ICdTcCcuDQoJLy9UaGVuIHdlIGludmVy dCAnU3AnIHRvICdTcHUnDQoJLy8tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0t LS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tDQoJZ3NsX21hdHJp eCAqIFNwID0gZ3NsX21hdHJpeF9hbGxvYyAobl9maWwsIG5fZmlsKTsNCgln c2xfbWF0cml4X3NldF96ZXJvIChTcCk7DQoJZm9yIChpID0gMDsgaSA8IG5f ZmlsOyBpKyspDQoJCWdzbF9tYXRyaXhfc2V0IChTcCwgaSwgaSwgZ3NsX3Zl Y3Rvcl9nZXQoUywgaSkpOwkvLyBWZWN0b3IgJ1MnIHRvIG1hdHJpeCAnU3An DQoJDQoJZ3NsX3Blcm11dGF0aW9uICogcCA9IGdzbF9wZXJtdXRhdGlvbl9h bGxvYyAobl9maWwpOw0KCWludCBzaWdudW07DQoJZ3NsX2xpbmFsZ19MVV9k ZWNvbXAgKFNwLCBwLCAmc2lnbnVtKTsJCQkJLy8gQ29tcHV0aW5nIHRoZSBM VSBkZWNvbXBvc2l0aW9uDQoJDQoJZ3NsX21hdHJpeCAqIFNJID0gZ3NsX21h dHJpeF9hbGxvYyAobl9maWwsIG5fZmlsKTsNCglnc2xfbGluYWxnX0xVX2lu dmVydCAoU3AsIHAsIFNJKTsJCQkJLy8gQ29tcHV0aW5nIHRoZSBpbnZlcnNl IHRocm91Z2ggTFUgZGVjb21wb3NpdGlvbg0KCWdzbF9wZXJtdXRhdGlvbl9m cmVlKHApOw0KCS8vZW5kIEludmVydGluZyBTLy8NCgkJDQoJCQ0KCQkNCgln c2xfbWF0cml4ICogVlQgPSBnc2xfbWF0cml4X2FsbG9jIChuX2ZpbCwgbl9m aWwpOw0KCWdzbF9tYXRyaXhfdHJhbnNwb3NlX21lbWNweSAoVlQsIFYpOwkJ CQkJLy8gVHJhbnBvc2Ugb2YgVg0KCQkNCgkJDQoJLy9USEUgUFNFVURPSU5W RVJTRS8vDQoJLy8tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0t LS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tDQoJLy9Db21wdXRhdGlvbiBvZiB0 aGUgcHNldWRvaW52ZXJzZSBvZiB0cmFucyhBKSBhcyBwaW52KEEpID0gVbdp bnYoUykudHJhbnMoVikgICB3aXRoIHRyYW5zKEEpID0gVS5TLnRyYW5zKFYp DQoJLy8tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0t LS0tLS0tLS0tLS0tLS0tLS0tDQoJZ3NsX21hdHJpeCAqIFNJcFZUID0gZ3Ns X21hdHJpeF9hbGxvYyAobl9maWwsIG5fZmlsKTsNCglnc2xfYmxhc19kZ2Vt bSAoQ2JsYXNOb1RyYW5zLCBDYmxhc05vVHJhbnMsCQkJCS8vIENhbGN1bGF0 aW5nICBpbnYoUykudHJhbnMoVikNCiAgICAgICAgICAgICAgICAJMS4wLCBT SSwgVlQsDQogICAgICAgICAgICAgICAgCTAuMCwgU0lwVlQpOw0KDQoJCQkN Cglnc2xfbWF0cml4ICogcGludiA9IGdzbF9tYXRyaXhfYWxsb2MgKG5fY29s LCBuX2ZpbCk7CS8vIENhbGN1bGF0aW5nICBVt2ludihTKS50cmFucyhWKQ0K CWdzbF9ibGFzX2RnZW1tIChDYmxhc05vVHJhbnMsIENibGFzTm9UcmFucywN CiAgICAgICAgICAgICAgICAJMS4wLCBVLCBTSXBWVCwNCiAgICAgICAgICAg ICAgICAJMC4wLCBwaW52KTsNCg0KCS8vZW5kIFRIRSBQU0VVRE9JTlZFUlNF Ly8NCgkNCiAgIAlzdGQ6OmNvdXQgPDwgInBpbnY6IiA8PCBzdGQ6OmVuZGw7 DQogICAJZm9yIChpID0gMDsgaSA8IG5fY29sOyBpKyspICANCiAgICAJCWZv ciAoaiA9IDA7IGogPCBuX2ZpbDsgaisrKQ0KICAgIAkJCXByaW50ZiAoIm0o JWQsJWQpID0gJWdcbiIsIGksIGosIA0KICAgIAkJCQlnc2xfbWF0cml4X2dl dCAocGludiwgaSwgaikpOw0KICAgCXN0ZDo6Y291dCA8PCAiXG4iIDw8IHN0 ZDo6ZW5kbDsNCgkNCg0KDQoJZ3NsX21hdHJpeF9mcmVlKFZUKTsNCglnc2xf bWF0cml4X2ZyZWUoU0kpOw0KCWdzbF9tYXRyaXhfZnJlZShTSXBWVCk7DQoJ Z3NsX21hdHJpeF9mcmVlKGdBX3QpOw0KCWdzbF9tYXRyaXhfZnJlZShVKTsN Cglnc2xfbWF0cml4X2ZyZWUoZ0EpOw0KCWdzbF9tYXRyaXhfZnJlZShWKTsN Cglnc2xfdmVjdG9yX2ZyZWUoUyk7DQoNCg0KDQoNCglyZXR1cm4gMDsNCn0= --Boundary-00=_tDz+HNGlP3bHTH2--