Hi everyone,
I am trying to find pseudoinverse of a matrix and as a first step I am trying to use svdcmp.c to find the SVD of the matrix. But I am unsure about how to call the function svdcmp() in my main since it uses **A (pointer to pointer) as one of its arguments. I dunno how to use them. Can anyone give me an example for finding SVD of a small 4x2 matrix using svdcmp() ?
thanks
shwe

Hey all,
I figured out the answer to my question and was also able to get the pseudoinverse using SVD.
Shwe

when i use Matlab to find Pseudo-Inverse Matrix of a matrix A(mxn), use svd()....A=U*S*Vt,.... , it's ok:

[m,n]=size(A);
[U,S,V]=svd(A);
r=rank(S);
SR=S(1:r,1:r);
SRc=[SR^-1 zeros(r,m-r);zeros(n-r,r) zeros(n-r,m-r)];
A_pseu=V*SRc*U.';

But now move to C++, i have a little problem with svd() in C. Cause , in matlab , U (mxm),S(mxn),V(nxn) ; in C++, U(mxn),S(nxn),V(nxn). The different is the sizes of those matrix U,S,V.

How could you find A_pseu in C ? :eek:

help me...

...in matlab , U (mxm),S(mxn),V(nxn) ; in C++, U(mxn),S(nxn),V(nxn).
How could you find A_pseu in C ?
Details are different; results are the same.

Example SVD problem http://en.wikipedia.org/wiki/Singular_value_decomposition. I'll find the pseudoinverse two ways.

I initialize A as follows, then use your matlab program. (Actually I use octave, since it's free.)


A = [ 1 0 0 0 2;0 0 3 0 0;0 0 0 0 0;0 4 0 0 0]
.
.
.
A_pseu =

0.20000 0.00000 0.00000 0.00000
0.00000 0.00000 0.00000 0.25000
0.00000 0.33333 0.00000 0.00000
0.00000 0.00000 0.00000 0.00000
0.40000 0.00000 0.00000 0.00000


Now, with Numerical Recipes SVD:

*****************Original Matrix****************
Matrix A:
1.000000 0.000000 0.000000 0.000000 2.000000
0.000000 0.000000 3.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 4.000000 0.000000 0.000000 0.000000

*******************After SVD********************
Matrix U:
0.000000 0.000000 1.000000 0.000000 0.000000
0.000000 -1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 1.000000
-1.000000 0.000000 0.000000 0.000000 0.000000

Diagonal of Matrix W:
4.000000 3.000000 2.236068 0.000000 0.000000

Matrix V:
0.000000 0.000000 0.447214 0.000000 0.894427
-1.000000 0.000000 0.000000 0.000000 0.000000
0.000000 -1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 1.000000 0.000000
0.000000 0.000000 0.894427 0.000000 -0.447214

************************************************
Matrix W
4.000000 0.000000 0.000000 0.000000 0.000000
0.000000 3.000000 0.000000 0.000000 0.000000
0.000000 0.000000 2.236068 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000

Matrix V-transpose:
0.000000 -1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -1.000000 0.000000 0.000000
0.447214 0.000000 0.000000 0.000000 0.894427
0.000000 0.000000 0.000000 1.000000 0.000000
0.894427 0.000000 0.000000 0.000000 -0.447214

*****Check product against original matrix******
Original matrix:
1.000000 0.000000 0.000000 0.000000 2.000000
0.000000 0.000000 3.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 4.000000 0.000000 0.000000 0.000000

Product of [U] * [W] * [V-Transpose]:
1.000000 0.000000 0.000000 0.000000 2.000000
0.000000 0.000000 3.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 4.000000 0.000000 0.000000 0.000000

************************************************
Matrix U-Transpose:
0.000000 0.000000 0.000000 -1.000000
0.000000 -1.000000 0.000000 0.000000
1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000

Matrix W+:
0.250000 0.000000 0.000000 0.000000 0.000000
0.000000 0.333333 0.000000 0.000000 0.000000
0.000000 0.000000 0.447214 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 0.000000

***********This is the Pseudo-Inverse***********
Product of [V] * [W+] * [U-Transpose]:
0.200000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.250000
0.000000 0.333333 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000
0.400000 0.000000 0.000000 0.000000


The first three items (U, Diagonal of W, and V) come from the Numerical Recipes SVD routine, everything else is from standard matrix manipulations. I showed the various intermediate results for my own amusement and edification.

Regards,

Dave

Footnote: The stuff I show in this post for the Numericsl Recipes solution came from the Numerical Recipes version 3 SVD class. Using the Numerical Recipes Version 2 C++ svdcmp() function had different (but equivalent) intermediate values but the same end results.

Bottom line: SVD decomposition itself is not unique (as mentioned in the text), but using it to obtain the pseudoinverse "works," the same way regardless of intermediate details.

Decomposition with different shapes of intermediate matrices still gives the same results. (When I say, "same results," from two different programs, I am glossing over the fact that implementation details may affect roundoff errors.)