Nathan - I've looked through the code a bit more and consulted my Scientific Computing book and concluded that the problem ultimately lies in the fact that the return values from the underlying SVD code can be very, very small. The PDL wrapper code attempts to make life easy for you by computing rotation matrices $u and $v, but it does not handle infinitesimal numbers properly. I am pretty sure that this sort of situation can be handled, if by nothing else than a check for the edge condition and some predefined behavior. However, I am not familiar with the underlying algorithm used and I am not sure precisely how to handle small values here.
However, I make two observations: (1) The code works for real data. :-) (2) The SVD code should be replaced with a more modern SVD algorithm. To illustrate my first point, note that the following (highly synthetic) situation pretty much works: --------------------%<-------------------- # Construct a weird A and print the results $A = pdl( [1, 2, 3], [1, 2 + 1e-15, 3]) print "A is:\n", $A - 2 # Calculate the svd ($u, $s_diag, $v) = svd($A); # Rebuild A from the svd $s = zeroes($v); $s->diagonal(0,1) .= $s_diag; $vt = transpose($v); $matcopy = $u x $s x $vt; # Print the rebuild matrix print $matcopy - 2 -------------------->%-------------------- A is: [ [ -1 0 1] [ -1 8.8817842e-16 1] ] matcopy is: [ [ -1 0 1] [ -1 8.8817842e-16 1] ] --------------------%<-------------------- Of course, SVD is SVD, and it should be solid as a rock. Although it's not as solid as you (or I) would like, I don't think we have any developers on the list who have both the knowledge and the inclination to implement a better method. After all, this is a volunteer effort. :-) If you find that this actually gives problems with real-life calculations, drop another note and you may be able to convince me or somebody else to re-examine the SVD calculations. David -- Sent via my carrier pigeon.
_______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
