Hi,

I'm a novice at NR programming, so forgive me if my problem is obvious to you! I've been trying to perform a SVD in Fortran90, with matrices as defined on the wikipedia page (http://en.wikipedia.org/wiki/Singular_value_decomposition) (I think this is different to the definition used by the svdcmp subroutine)

I've followed the method in the Riley, Hobson and Bence textbook (which finds the eigenvalues and eigenvectors of the two matrices A*(AT) and (AT)*A) and used the jacobi and eigsrt subroutines in the relevant places. My problem is that both my U and VT matrices come out with the wrong elements, and so when I recombine U*S*VT, I don't get the original matrix. At the end of this post are the answers from the textbook.

Would anyone be able to explain how I fix up U and V? Or, could they check if using svdcmp gives the right answer, as I've not been able to make it do so!

Many thanks

Ravi

A:
2 2 2 2
1.7 0.1 -1.7 -0.1
0.6 1.8 -0.6 -1.8

SVD should give...

U:
1.0 0.0 0.0
0.0 0.6 -0.8
0.0 0.8 0.6

S:

4 0 0 0
0 3 0 0
0 0 2 0

VT:

0.5 0.5 0.5 0.5
0.5 0.5 -0.5 -0.5
-0.5 0.5 0.5 -0.5
0.5 -0.5 0.5 -0.5

check if using svdcmp gives the right answer



How did you fill the matrix and call svdcmp? Here's what I got from GNU gfortran on my Linux system, using Numerical Recipes version 2.1 svdcmp.f90 and friends:


Original Matrix U:
2.000000 2.000000 2.000000 2.000000
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600000 -1.800000

After svdcmp(u(1:m,1:n), s(1:n), v(1:n,1:n) with m = 3 and n = 4

Matrix U:
-1.000000 0.000000 0.000000 0.000000
0.000000 0.600000 0.000000 0.800000
0.000000 0.800000 0.000000 -0.600000

Diagonal of Matrix S:
4.000000 3.000000 0.000000 2.000000

Matrix V:
-0.500000 0.500000 -0.500000 0.500000
-0.500000 0.500000 0.500000 -0.500000
-0.500000 -0.500000 -0.500000 -0.500000
-0.500000 -0.500000 0.500000 0.500000

V-Transpose:
-0.500000 -0.500000 -0.500000 -0.500000
0.500000 0.500000 -0.500000 -0.500000
-0.500000 0.500000 -0.500000 0.500000
0.500000 -0.500000 -0.500000 0.500000

Product U * S * [V-Transpose]
(Compare with original)
2.000001 2.000001 2.000000 2.000001
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600001 -1.800000



Regards,

Dave

I filled the matrix in within the program, since I just want to test it yet


A(1,1) = 2; A(1,2) = 2; A(1,3) = 3; A(1,4) = 2
A(2,1) = 1.7; A(2,2) = 0.1; A(2,3) = -1.7; A(2,4) = -0.1
A(3,1) = 0.6; A(3,2) = 1.8; A(3,3) = -0.6; A(3,4) = -1.8


To call the subroutine, I just used call svdcmp(A,w,v), and I used the svdcmp_dp subroutine. The output I get is this:


U
-0.9909 0.1202 0.0000 0.0601
0.1188 0.5755 0.0000 0.8091
0.0627 0.8089 0.0000 -0.5846
W
4.60888 2.96958 6.942025E-25 1.98490
Vt
-0.3780 -0.4030 -0.6970 -0.4571
0.5738 0.5906 -0.3715 -0.4288
-0.4417 0.5522 -0.4417 0.5522
0.5768 -0.4288 -0.4255 0.5499
Recomb A
2.0000 2.0000 3.0000 2.0000
1.7000 0.1000 -1.7000 -0.1000
0.6000 1.8000 -0.6000 -1.8000


As you can see, the wrong intermediate matrices come out, but they combine to give the right answer! I'm not entirely sure what I'm doing wrong...:confused:

I've put the svdcmp.f90 in a module but when I try to build the file, I get the following error message

Access violation:
The instruction at address 00453ba5 attempted to read from location 00000014
00453b83 describe_module_routine_to_linker(<ptr>structÄscoped_entity)#1D [+0022]

004582be read_module_entity(enumÄlogical)#1D [+12b3]

00459a96 process_binary_module(int,enumÄlogical)#1D [+034b]

0045a282 process_use_stmt(<ptr>char,<ref>int) [+0264]

0040a953 parse_declaration_statement(<ptr>char,int,int,<ref>int) [+25e2]

0041202b handle_token(<ptr>char,int,int,int,int,<ref>int) [+0df5]

0041a8f8 get_implicit_interfaces(void)#B [+0748]

0041bbd3 process_contained_interfaces(void) [+00cc]

eax=00000000 ebx=00000000 ecx=0383f1ec
edx=00000000 esi=00000003 edi=0383f288
ebp=0383f12c esp=0383ed08 IOPL=0
ds=0023 es=0023 fs=003b
gs=0000 cs=001b ss=0023
flgs=00210202 [NC OP NZ SN DN NV]
0360/6820 TSTK=5 [ ]

00453ba5 mov edi,[eax+0x14]

00453ba8 mov ecx,0xffffffff

A(1,1) = 2; A(1,2) = 2; A(1,3) = 3;<---Different from previous A(1,4) = 2
A(2,1) = 1.7; A(2,2) = 0.1; A(2,3) = -1.7; A(2,4) = -0.1
A(3,1) = 0.6; A(3,2) = 1.8; A(3,3) = -0.6; A(3,4) = -1.8


As you can see, the wrong intermediate matrices come out No, actually I can't see that. Your previous example (and the one for which I posted my solution) had A(1,3) = 2.

but they combine to give the right answer! Well, with significantly different intermediate values, that would be a head-scratcher, for sure! However, see Footnote.

Regards,

Dave

Footnote:
"The universe is not only stranger than we imagine; it is stranger than we can imagine."
---Arthur C. Clarke (1917-2008)

I've put the svdcmp.f90 in a module but when I try to build the file, I get the following error message
How, exactly, did you change svdcmp.f90 to make "put it in" a module?
How, exactly, did you compile that file?
How, exactly, did you compile and link your test program? (Needs stuff from pythag and nrutil, right?)

Regards,

Dave

I'm such a goon! The reason U, S, V looked different is because of that wrong element! I didn't even notice it until you mentioned it. Once I'd corrected it, everything worked fine!

Thanks muchly for your help...I was starting to despair, but now I can get on with my project!

hi
i have quastion about svd ,but i do not how can i ask my quastion???
please help me!
regards

...help...
If your question is about how to use or how to test some functions defined in Numerical Recipes, then just tell us what you don't understand.

If you have a question about something in this thread, then just ask. If it is about C or C++ rather than Fortran, you can tell us what language (what compiler) you are using.

If your question is not directly related to the topics already covered in this thread, then you can start another thread.

Ask specific questions.

Regards,

Dave

hi
i need w matrix organize big to small.but the svdcmp fortran subroutine sometimes give w in irregular way.
can you tell me why this problem occur?(because w must be organized bog to small w1>w2>w3....)
and can you tell me how can i solve this problem?
please give help in emergency!!
thanks

hi
i need w matrix organize big to small...

Suppose you have the following, from xsvdcmp.f90, as I showed in the previous example in this thread:


Decomposition Matrices:
Matrix U
-1.000000 0.000000 0.000000 0.000000
0.000000 0.600000 0.000000 0.800000
0.000000 0.800000 0.000000 -0.600000
Diagonal of Matrix W
4.000000 3.000000 0.000000 2.000000
Matrix V-Transpose
-0.500000 -0.500000 -0.500000 -0.500000
0.500000 0.500000 -0.500000 -0.500000
-0.500000 0.500000 -0.500000 0.500000
0.500000 -0.500000 -0.500000 0.500000
Check product against original matrix:
Original Matrix:
2.000000 2.000000 2.000000 2.000000
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600000 -1.800000
Product U*W*(V-Transpose):
2.000001 2.000001 2.000000 2.000001
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600001 -1.800000

Now, if you want the Singular Values (the elements if vector W) to be sorted, then here's the drill:

If you want to exchange the third and fourth elements of W, then, since these represent the diagonal elements of the diagonal matrix in the middle of the product [U] * [W] * [V-transpose] you should exchange the third and fourth columns of U and exchange the third and fourth rows of V-Transpose (the third and fourth columns of V that svdcmp gives).

Sometimes people take the trouble to get the maximum number of non-negative values in U and V by inverting signs of columns of U and corresponding rows of V-transpose as needed.

Final results might be something like

*************************************
Sorting the Singular Values
Swapping col 4 and 3 of U and V
*************************************
After sorting W and inverting signs
Matrix U
1.000000 0.000000 0.000000 0.000000
0.000000 0.600000 -0.800000 0.000000
0.000000 0.800000 0.600000 0.000000
Diagonal of Matrix W
4.000000 3.000000 2.000000 0.000000
Matrix V-Transpose
0.500000 0.500000 0.500000 0.500000
0.500000 0.500000 -0.500000 -0.500000
-0.500000 0.500000 0.500000 -0.500000
0.500000 -0.500000 0.500000 -0.500000
Check product against original matrix:
Original Matrix:
2.000000 2.000000 2.000000 2.000000
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600000 -1.800000
Product U*W*(V-Transpose):
2.000001 2.000001 2.000000 2.000001
1.700000 0.100000 -1.700000 -0.100000
0.600000 1.800000 -0.600001 -1.800000


Regards,

Dave

thanks alot.
but i think in solution u(1,1)=-1.000000
thanks

hi
in solution u matrix and v matrix in my run differ in sign at some columns. but product of u*w*v(transpose) is valid.
can you tell me why in spite of this diference product of u*w*v(transpose) is valid?
and i run svdcmp as you but why this difference appeare?
mybe my subroutine is not correct ????!!!!!!

in solution u matrix and v matrix in my run differ in sign at some columns. but product of u*w*v(transpose) is valid. As a specific example: If you change the signs on all elements of column number 3 of matrix U and column number 3 of matrix V (that is, row number 3 of matrix V-Transpose), the product of the three is the same. It's fundamental.
...why this difference...The decomposition is not unique. Different implementations and different versions may come up with different solutions. Furthermore, depending on the compiler and roundoff error, small values (absolute value near zero) may come out as zero or "minus zero", giving a different answer when the decision is made to change signs of a particular column. If you print with E format instead of F, you might see some differences in the results.
...mybe my subroutine is not correct...If the Singular Values are the same and the product U times W times V-Transpose gives the original matrix, I would say that your subroutine is probably correct. Try different sizes and shapes of examples, and try example matrices with various condition numbers. If you get the original matrix back each time (or close to it, taking roundoff errors into account), it's probably OK.

I mean: In general you can't prove a program like this is correct by testing, but testing with lots of different kinds of matrices (and knowing how to interpret the results) can usually give you some confidence in it.

Regards,

Dave

hi
before i use svdcmp i use from imsl library for svd .i use (lsvrr) but when LDA is grater than around (500) solution was wrong but for lda later than mentined solution was vali .can you tell me why?????

hi
in my running of svdcmp the data and result are:
original matrix
2 -------- 2 -------- 2 -------- 2
1.7 -------- .1 -------- -1.7 -------- -.1
.6 -------- 1.8 -------- -.6 -------- -1.8

w

4 -------- 3 -------- 2 -------- 4.67E-29

u matrix

-1 -------- -1.73E-17 -------- 3.8E-16 -------- 3.86E-32
-2.21E-17 -------- .6 -------- -.8 -------- -3.81E-30
2.41E-16 -------- .8 -------- .6 -------- -5.09E-30

v matrix

-.5 -------- .5 -------- -.5 -------- .5
-.5-------- .5 -------- .5 -------- -.5
-.5 -------- -.5 -------- .5 -------- .5
-.5 -------- -.5 -------- -.5 -------- -.5


original matrix(u*w*v(trans))

2 -------- 2 -------- 2 -------- 2
1.7 -------- .1-------- -1.7 -------- -.1
.6 -------- 1.8 -------- -.6-------- -1.8

you can see first column in u and v are difer in sign with your solution.
my question is if we are run same program why results can have diferences?

...f we are run same program ...

1. My previous output was run with single precision variables. Your output is apparently with double precision variables.

2. Who says that we are running the same program? I added the column-swapping and sign flipping to my test program after getting the results from svdcmp. I did mention that results (decision about which columns would undergo a sign change) may vary with compiler and with roundoff errors inherent in any floating point calculations.

For example, here is your U matrix output, where I have reformatted it to look like something that I might want to read:

Your U matrix:
-1.000000E+00 -1.730000E-17 3.800000E-16 3.860000E-32
-2.210000E-17 6.000000E-01 -8.000000E-01 -3.810000E-30
2.410000E-16 8.000000E-01 6.000000E-01 -5.090000E-30


Here's the output that I got from Numerical Recipes Version 2.10 svdcmp.f90 when I ran it with double precision variables, using GNU gfortran 4.1.2, and when I used the same E format:

Matrix U
-1.000000E+00 -6.893234E-16 2.688982E-32 3.456942E-16
-3.996947E-16 6.000000E-01 -2.656239E-30 -8.000000E-01
-2.027172E-16 8.000000E-01 -3.541653E-30 6.000000E-01


Here is matrix U from my test program after my test program swapped columns three and four and flipped the signs of the first column. (It also swapped vector W elements 3 and 4 and V matrix columns 3 and 4 and flipped the signs of V matrix column 1. I only show the U matrix here.)


Matrix U
1.000000E+00 -6.893234E-16 3.456942E-16 -2.688982E-32
3.996947E-16 6.000000E-01 -8.000000E-01 2.656239E-30
2.027172E-16 8.000000E-01 6.000000E-01 3.541653E-30



I don't see any significant discrepancies. If you are demanding exact (bit-for-bit) results in any floating point calculations, and you don't accept inevitability of roundoff errors that can affect certain decisions, you will surely be disappointed. Note, particularly that for your results and for mine too, the product of [U] * [W] * [V-Transpose] is equal to the original matrix (plus or minus a relatively small roundoff error).

With double precision variables, the final product was the following (whether I swapped columns or flipped signs or not):

Product U*W*(V-Transpose):
2.000000E+00 2.000000E+00 2.000000E+00 2.000000E+00
1.700000E+00 1.000000E-01 -1.700000E+00 -1.000000E-01
6.000000E-01 1.800000E+00 -6.000000E-01 -1.800000E+00


So, it can be seen that this agrees with the original matrix to seven significant decimal digits.

Bottom line: (These are all rhetorical questions. Don't tell me the answers. Ask yourself the questions and proceed accordingly.)

Why are you performing the decomposition?

What is your application program that uses the SVD results?

Did your SVD results give you something that "worked" in the sense that the final output met your expectations?

Is the outcome of your application program affected by the small differences in roundoff error between your SVD results and my test program's SVD results? (If it is, then what?)

Regards,

Dave

"The purpose of computing is insight not numbers"
---Richard W. Hamming

HI
thanks for previous answers.
if i use double percision what is error in array evaluting of u or v matrix?
i need know what is error in evaluting u or v?for example(e=1E-008)
thanks

hi
i want to solve linear equation(Ax=B).
this equation is not solve with ordinary method.
first I thought A matrix is singular but when I used svd for A
matrix I recognized that this is well-conditioned (rank of A=30).
and when I saw U and V matrix I saw that some array of them
are for example(1.5E-020,etc.).so I think problem origin from
u and v? so I need error in evaluting U and V?
can you help me in emergency?
regard
amin

hi
I do not know why you do not ask my previous question.(I need the answers in emergency).please ask following question:
in evaluating columns of U or V we know the columns divide by norm of the column.I want to know the norm of culumn befor the norm become (1.0).please give me answer???(I use svdcmp subroutine).....please answer me............(I need.........)
regards
amin

...I need...
I didn't mean to be rude by not answering your most recent questions or meeting your demands. If I thought I could help further, I would try.

Sorry.

Maybe someone else...

Regards,

Dave

Hi all,

plaise can you help me to compute the complex SVD in Matlab

thanks in advance