Skip to content

SVD routines sometimes return -0 as a singular value #324

Open
@t-mitchell

Description

@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

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions