On Tue, Jun 22, 2010 at 2:06 PM, David Mertens <[email protected]>wrote:
> 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. > Also, feel free to file this as a bug. I would give it a title such as "SVD does not handle severely pathological cases" or some such. David -- Sent via my carrier pigeon.
_______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
