Svd
Has anyone else used the svdcmp.c routine and seen abysmal speed performance? I am using matrices on the order of 500x500 (double percision) and the NRC routine is taking > 10 seconds (compiled with Borland C++ Builder v6).
The exact same matrices in Mathematica (v4.1) using the built-in function SingularValues[] take < 2 seconds.
One reason perhaps has to do with sparsity; my matrices are 90+% sparse, and Mathematica might be handling this much more efficiently.
Still, I am floored by the difference in performance between compiled NRC code and Mathematica.
My app runs in realtime, and I cannot afford 10 seconds.
Any suggestions?
Thanks in advance.
John W. Szuhay
Houston, Texas
Yikes...I just discovered that the NRC svd routine is taking > 100 iterations to solve my problem!
I am normalizing each column before I send the matrix to svdcmp, and wonder now if Mathematica is running so much faster because of better scaling.
What factors affect the speed of convergence in svdcmp?
I do not know the details of the algorithm, but I suspect the technique may be prone to numerical instability for some cases...
BTW, the ratio of largest/smallest singular value is ~1.7/0.0023.
Also, Mathematica and svdcmp give the same singular values.
Any SVD wizards araound here?
Thanks,
John
Methinks the best way for you to learn about SVD is to either read the original article by Golub and Reinsch (Numer. Math. 14, 403-420 (1970)), whose included ALGOL program is the basis for the svdcmp() routine, or the book "Matrix Computations" by Golub and Van Loan.
You might also be interested in T.F. Chan's SVD algorithm (ACM algorithm 581 (http://www.netlib.org/toms/581) ).
Suffice it to say that Mathematica is specially programmed to handle sparse matrices efficiently, hence the time difference you have said.
I have not yet done any research on more efficient methods for computing the SVD of a sparse matrix, but you can try looking on the Web.
Hope this helps.
Jan M. (^_^)
P.S. The SVD algorithm is (theoretically) stable. See the above references for an explanation.
I decided to dig a little deeper and discovered:
-> the original svdcmp.c routine in NRC has IMHO a very poor convergence strategy; when I relaxed the convergence criteria I got a signiiicant speedup without loss of accuracy (for my problem).
-> the sparse svd routines available on the net such as svdpackc, svdlibc. etc. do not perform any faster than svdcmp.c.
Unfortunately I am still stuck. I am trying to invert a 300x300 and want to use svd since this array can become very close to singular.
BTW, my matrix is 0.90% dense. :eek:
John
Exactly how did you "relax" the convergence criteria?
Jan M. ('_')