I'm converting an expectation maximization algorithm for mixtures of PCA analyzers from Matlab to C++ and everything was going fine until I hit the SVD routine. When the matrices are small, NR and Matlab spit up the same results; however, when the matrices get large enough, the columns of the U matrix corresponding to the zero eigen values are very small when using the NR algorithm but not for the Matlab algorithm. I have an orthogonality check which returns true for both algorithms but the NR method is only succeeding because the column elements are all very close to zero and not because it is producing truly orthogonal columns. You may be wondering why I care about orthogonality in the NULL space but the maximization step depends heavily on this. If any of you care to see this for yourself, just choose any square matrix whose size is greater than 20 but is not full rank. Then take the SVD in matlab and the SVD in NR and the columns of U corresponding to the NULL space will NOT match. Can anybody tell me what's happening?

Original matrix:
-0.0109 0.0342 -0 -0.0135 -0.0045 -2.9321 -0.1268 1.3343 1.7193
1.0095 1.4301 0 0.0443 0.2424 -1.86 -0.1874 0.4193 -1.0982
-0.9743 -0.6319 -0 0.2606 0.0073 -1.473 0.9752 1.0572 0.7789
-0.5112 -0.6011 -0 2.5637 -0.0427 0.5127 -0.301 -0.3034 -1.317
-0.3558 0.8159 -0 -0.5747 0.3286 2.8233 -1.0772 -0.8434 -1.1167
0.7852 1.0328 0 -0.3388 -2.1324 -0.9182 0.6886 0.2664 0.6165
0.395 1.6574 0 0.2399 -2.0996 -1.0576 0.2325 0.3509 0.2815
0.9565 1.6878 0 -0.0367 -0.163 -1.1683 -0.879 -0.0839 -0.3134
-0.3583 -1.2252 0 0.8907 0.0368 -0.7362 -0.1266 1.2813 0.2375
0.8135 0.7549 -0 0.4182 -0.2464 -1.5103 0.614 -1.4011 0.5573
0.7765 1.1075 0 0.2711 1.0047 -0.8255 0.0343 -1.52 -0.8487
0.2276 0.3095 0 0.9551 0.7178 1.4594 -0.2679 -1.6561 -1.7454
-0.3967 0.0867 -0 -1.2165 0.7947 -2.0118 0.9085 -0.1303 1.9653
-0.8158 0.2208 -0 0.3796 0.2872 2.9451 -0.9654 -0.8962 -1.1553
0.1617 0.3088 0 -2.4751 1.615 -0.5028 -0.119 0.2142 0.797
-0.4721 -0.9035 -0 1.0502 1.0622 1.4735 0.2446 -1.7523 -0.7026
-1.0334 -1.8056 -0 -0.211 0.3782 1.3945 0.4733 0.4685 0.3355
-1.0857 -1.3437 -0 -0.5717 0.0866 1.3886 0.4473 0.0253 1.0533
1.0307 1.5281 0 -0.0367 -0.0313 -1.6544 -0.2569 -0.3498 -0.2297
-1.3613 -0.478 -0 1.58 0.3853 0.969 -1.1876 0.5244 -0.4318


NR's results:
Matrix U
0.370428 0.164397 -0.144067 0.144678 -0.209874 0.325576 0.265529 -0.103779 0.00122625
0.191291 -0.293723 -0.0514942 0.13395 -0.307787 -0.306485 -0.363504 -0.315551 -0.305766
0.164606 0.250351 -0.1878 0.126281 0.0283952 0.03901 -0.666419 -0.26169 0.205138
-0.196207 -0.095839 -0.483822 0.207678 0.0210748 0.0326778 -0.0740954 0.350909 -0.113084
-0.324222 -0.166328 0.276962 -0.232425 -0.182776 0.231088 -0.117418 -0.227012 0.467981
0.216527 -0.108012 -0.0772437 -0.439896 0.3001 -0.0466633 -0.0257833 -0.0417638 -0.433948
0.20355 -0.195981 -0.182655 -0.420408 0.143881 0.243618 -0.288861 0.4557 0.169896
0.147513 -0.302801 0.0234165 -0.0244385 -0.240092 0.185758 0.203154 -0.0347863 -0.175546
0.0476944 0.205516 -0.304748 0.109205 -0.198343 -0.156546 0.206152 0.141696 0.0814414
0.168409 -0.19085 -0.00521466 0.19344 0.44437 0.221006 0.218386 -0.376316 0.119145
0.0273787 -0.320746 0.115534 0.369756 0.0988074 0.0291716 -0.130121 0.335944 0.197659
-0.27053 -0.280931 0.00805484 0.207962 0.116494 -0.0863898 -0.0968139 -0.0140165 -0.052954
0.280503 0.199581 0.252034 0.252798 0.180052 0.391962 -0.242678 0.0888883 -0.200142
-0.383338 -0.104676 0.0833831 -0.14038 -0.0951167 0.314055 -0.130931 -0.0696804 -0.441619
0.120712 0.112933 0.555187 0.137448 -0.286169 -0.0320187 -0.0591723 0.367293 -0.15616
-0.280459 -0.0240895 -0.0033157 0.350816 0.32531 0.0703665 0.0321855 0.00411801 -0.177501
-0.167363 0.346348 0.0159304 -0.00411212 0.0611179 -0.185903 -0.0856765 -0.0593264 -0.125899
-0.127106 0.334234 0.121681 -0.093988 0.196823 0.159042 0.026107 -0.0020747 -0.0791541
0.194705 -0.288341 0.0281673 0.0892633 -0.072195 0.0708231 0.0384902 -0.0781372 0.00751603
-0.199738 0.0775624 -0.299621 0.0739627 -0.34806 0.493221 -0.04493 -0.0619176 -0.09374
Matrix S
8.77693 0 0 0 0 0 0 0 0
0 6.24593 0 0 0 0 0 0 0
0 0 4.78824 0 0 0 0 0 0
0 0 0 3.94075 0 0 0 0 0
0 0 0 0 3.154 0 0 0 0
0 0 0 0 0 1.89976 0 0 0
0 0 0 0 0 0 1.21861 0 0
0 0 0 0 0 0 0 5.5773e-05 0
0 0 0 0 0 0 0 0 5.10836e-16

Matlab's results:
U =

0.3704 0.1644 0.1441 -0.1447 -0.2099 0.3256 0.2655 0.1038 0.2743
0.1913 -0.2937 0.0515 -0.1339 -0.3078 -0.3065 -0.3635 0.3156 0.5169
0.1646 0.2504 0.1878 -0.1263 0.0284 0.0390 -0.6664 0.2617 -0.5458
-0.1962 -0.0958 0.4838 -0.2077 0.0211 0.0327 -0.0741 -0.3509 0.0575
-0.3242 -0.1663 -0.2770 0.2324 -0.1828 0.2311 -0.1174 0.2270 -0.0047
0.2165 -0.1080 0.0772 0.4399 0.3001 -0.0467 -0.0258 0.0418 0.0431
0.2036 -0.1960 0.1827 0.4204 0.1439 0.2436 -0.2889 -0.4557 0.1259
0.1475 -0.3028 -0.0234 0.0244 -0.2401 0.1858 0.2032 0.0348 -0.3513
0.0477 0.2055 0.3047 -0.1092 -0.1983 -0.1565 0.2062 -0.1417 -0.1381
0.1684 -0.1908 0.0052 -0.1934 0.4444 0.2210 0.2184 0.3763 -0.0813
0.0274 -0.3207 -0.1155 -0.3698 0.0988 0.0292 -0.1301 -0.3359 -0.0372
-0.2705 -0.2809 -0.0081 -0.2080 0.1165 -0.0864 -0.0968 0.0140 0.0408
0.2805 0.1996 -0.2520 -0.2528 0.1801 0.3920 -0.2427 -0.0889 0.2884
-0.3833 -0.1047 -0.0834 0.1404 -0.0951 0.3141 -0.1309 0.0697 0.0460
0.1207 0.1129 -0.5552 -0.1374 -0.2862 -0.0320 -0.0592 -0.3673 -0.1447
-0.2805 -0.0241 0.0033 -0.3508 0.3253 0.0704 0.0322 -0.0041 0.0441
-0.1674 0.3463 -0.0159 0.0041 0.0611 -0.1859 -0.0857 0.0593 0.2280
-0.1271 0.3342 -0.1217 0.0940 0.1968 0.1590 0.0261 0.0021 0.0859
0.1947 -0.2883 -0.0282 -0.0893 -0.0722 0.0708 0.0385 0.0781 -0.1242
-0.1997 0.0776 0.2996 -0.0740 -0.3481 0.4932 -0.0449 0.0619 0.0597


S =

8.7769 0 0 0 0 0 0 0 0
0 6.2459 0 0 0 0 0 0 0
0 0 4.7882 0 0 0 0 0 0
0 0 0 3.9408 0 0 0 0 0
0 0 0 0 3.1540 0 0 0 0
0 0 0 0 0 1.8998 0 0 0
0 0 0 0 0 0 1.2186 0 0
0 0 0 0 0 0 0 0.0001 0
0 0 0 0 0 0 0 0 0.0000

It seems support from the authors in these forums is non-existent and my struggle with their poorly implemented SVD routine has been frustrating to say the least. I chose to use the NR code because I didn't want to learn all the details of the SVD algorithm. Instead, the exact opposite has happened and I've learned that the NR SVD routine performs very poorly with REAL world HIGH dimensional data when compared to Matlab which uses LAPACK. I'm very discouraged and I recommend anyone who needs a real world SVD to use LAPACK. Yes it is platform dependent and a pain in the butt to link to externally but you won't waste 1 month of valuable research time figuring out why your "black-box" SVD routine is spitting out non-orthogonal singular vectors. I only hope that there is SOME valuable content in the rest of the NR book/code.

Hi,

The example you gave of a matrix that performs poorly with svdcmp.cpp seems fine to me. The columns of U are orthogonal, and it's got nothing to do with having small elements. In fact, the smallest element is of order .001, and the number
of such small elements in the Matlab output is roughly the same.

The reason that Matlab gives different output is because the SVD is not unique, as discussed in the text.

It is true that if you have *very* badly scaled data in your input matrix, then the algorithm in Lapack can often do a better job than the algorithm in NR. However, the example you gave is not of this kind.

By the way, your comment about lack of support from the authors is a little unfair. For what you spent on the book and the software, the best you can expect is that they fix actual bugs. Commercial software that includes consulting is extremely expensive, because someone has to pay for people's time.

i hope you havent forgotten to Zero out the near zero elements in S so they dont have the disturbing effects on the S^-1 in the solving of problems

Sorry for the confusion, the example I provided above was only to demonstrate the manner in which the Matlab result and NR result differ. I wasn't trying to demonstrate a non-orthogonal case. In order for me to demonstrate that, I would have to use much higher dimonsional data and nobody wants me filling up this forum with numbers. If you would like, I will email you the data in a tar file so that you can see what I mean. The data is not badly scaled from what I can tell. I used the example above because I found it some what strange that the columns would match almost perfectly except for the last one, the one corresponding to the lowest singular value.

Also, I do not think my previous comment was unfair. I was probably a little too harsh, but not unfair. I spent about $70 on the book, $35 on the example book, and $120 on the software totalling $225. I have received better support from Microsoft's crappy OS and it cost me only $89. I have received better support from many other software packages costing only $50. For $225, I think I deserve some kind of response, even if that response is "leave me alone, I don't know."

Finally, it is important to remember that I didn't want to know the SVD details, I simply wanted an algorithm that gave me the "best" possible result. NR, in my opinion, does not do that.

Once again, thank you for your reply.

thanks, but I have done this in my actual application.

hi
if you like you can send me the matrix
i am at moment implementing the svd new :O)))
cos of i have some probs with it too
first its still the same as the one golumb and reinsch described.

otherwise would say you should try to solve problems with it and see if the different implementiaons behave really different in solving.

i have discovered that the "normal" NR SVD behaves very well in practical applications
in speed and accuracy

just thtere are some cases i think there is not enought flexibility for the problem.
the problem is mstly in small rotations that are not possible to perform cos of a large and a very small eigenvalue.

most of the probs have occured cos i decided to to

At*A= V*St*S*Vt

:O)