[Octave-bug-tracker] [bug #55710] svd may sometimes return -0 for exactly zero singular values

[Tim Mitchell]

URL:
  <https://savannah.gnu.org/bugs/?55710>

                 Summary: svd may sometimes return -0 for exactly zero
singular values
                 Project: GNU Octave
            Submitted by: tmitchell
            Submitted on: Thu 14 Feb 2019 08:33:10 AM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Tim Mitchell
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 4.0.0
        Operating System: Any

    _______________________________________________________

Details:

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 (negative 0).  While
allowable in IEEE floating point, returning -0 breaks the assumption that all
singular values should be nonnegative.  It would much better to ensure the
sign bits of all returned singular values are positive.  Otherwise, for
example, if the reciprocals are needed, one will get -inf for -0, which is not
at all intuitive.

This -0 behavior occurs in both Octave and MATLAB.  It appears to be coming
from the LAPACK SVD routines (but this isn't confirmed yet).  See [url
https://github.com/Reference-LAPACK/lapack/issues/324 this issue] for more
details.  

Here is a demo (which I've also attached since it includes comments):


    n = 10;
    A = randn(n);
    A = triu(A);
    
    A(1) = 0;
    fprintf('\nA = randn(n); A = triu(A); A(1) = 0\n')
    fprintf('Reciprocal of smallest singular value:\n');
    s = svd(A);
    fprintf('s = svd(A):       %g\n',s(end));
    [U,S,V] = svd(A);
    fprintf('[U,S,V] = svd(A): %g\n\n',S(end));

    A(1) = -0;
    fprintf('Same matrix but now A(1) = -0\n');
    fprintf('Reciprocal of smallest singular value:\n');
    s = svd(A);
    fprintf('s = svd(A):       %g\n',s(end));
    [U,S,V] = svd(A);
    fprintf('[U,S,V] = svd(A): %g\n\n',S(end));


Note that the -0's may occur in either case: A(1) = 0 or A(1) -0.  If n is set
>= 26, GESDD will be called instead of GESVD and so the minimum singular value
will not be exactly zero; hence for n >= 26, the demo is inconclusive.

Even 1/svd(-0) returns -inf (Octave 4.0.0 and R2017a).  However, on R2018, I
get +inf for 1/svd(-0).





    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: Thu 14 Feb 2019 08:33:10 AM UTC  Name: svdzero.m  Size: 2KiB   By:
tmitchell

<http://savannah.gnu.org/bugs/download.php?file_id=46251>

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?55710>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/