SVD routines sometimes return -0 as a singular value

[t-mitchell]

If exactly zero is returned as a singular value (unlikely but still happens sometimes), it may be that it is returned as +0 or -0. While allowed in IEEE floating point, returning -0 breaks the assumption that all singular values should be nonnegative. For example, if their reciprocals are needed, one would get -inf instead of +inf for this -0 singular value, which is not at all intuitive. It would be much better if the SVD routines ensured all singular values have positive sign bits, including 0. Fortunately, this is easy to do. Adding the rather necessary comment to explain why the absolute value of singular values is taken is probably more work. ;-)

I’ve seen this -0 issue in both MATLAB and Octave so I assume this is coming from LAPACK (but this should be confirmed directly). My attached MATLAB demo (rename the extension from .txt to .m) makes a random upper triangular matrix A of size n and sets A(1) = 0 and then calls svd in various ways. It repeats the tests with A(1) = -0. In either case, svd may sometimes return -0. In the demo, this happens for GESVD when singular vectors are also requested. On R2017b, even 1/svd(-0) returns -inf (but on R2018a it returns +inf). I haven't seen GESDD return exactly zero as a singular value so I can't confirm whether or not GESDD also has this -0 issue.

For an application where +0 versus -0 can be an issue, consider pseudospectra: the set of all eigenvalues of A + E, where norm(E) <= epsilon. A point z is inside the pseudospectrum if and only if the reciprocal of the smallest singular of zI - A is >= 1/epsilon. Clearly eigenvalues of A should be in the pseudospectrum. But if z is an eigenvalue of A, then the minimum singular value should be zero and its reciprocal should be +inf. But if -0 is returned, one will get -inf as the reciprocal and the test will say that eigenvalue z of A is actually not in the pseudospectrum.

svdzero.txt

[langou]

Hi Tim, this is a good point. I think it is an easy enough fix and makes total sense to me. Thanks for contributing! Cheers, Julien.

[t-mitchell]

As a more cautious alternative to taking the absolute value before returning the singular values, how about just explicitly setting any exactly zero singular value to +0 before returning? The argument is if there's ever a bug down the road that causes nonsensical singular values (e.g. negative), some one might notice it. With the absolute value, we might never know.