bugGNU Octave - Bugs: bug #55710, svd may sometimes return -0 for...

 
 

bug #55710: svd may sometimes return -0 for exactly zero singular values

Submitter:  Tim Mitchell <tmitchell>
Submitted:  Thu 14 Feb 2019 08:33:10 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Tim Mitchell Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 15 Feb 2019 01:02:28 AM UTC, comment #7: 

@Rik,

Sorry to cause you extra work, but is better solution.
And, may avoid issues later.

Michael Godfrey <godfrey>
Group Member
Fri 15 Feb 2019 12:49:41 AM UTC, comment #6: 

No, it's a fine idea.  I just coded up quickly using abs() because it is a built-in method on top of Octave Array objects in C++.  It turns out svd is implemented in C++ so I can't use the Octave language [S(S==0) = 0] to do this, but I can code a for loop for it.  I checked in this change here (https://hg.savannah.gnu.org/hgweb/octave/rev/75d79c39ac92).

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Thu 14 Feb 2019 11:07:17 PM UTC, comment #5: 

@Rik,

Was Tim's S(S==0)=0 not a good idea?

Michael Godfrey <godfrey>
Group Member
Thu 14 Feb 2019 06:44:09 PM UTC, comment #4: 

@Rik:

1) I agree that this -0 is likely to be an upstream issue with LAPACK.  Fortunately, it's easy to fix so I imagine it won't be too long from them to handle it on their end.

2) In the meantime, having Octave implement a temporary workaround would be ideal.  It's debatable which is better but I think S(S==0)=0 is better than abs(S).  I doubt there's a difference speedwise.  But if there's ever a bug where the SVD routines return a negative singular value, with S(S==0)=0, that negative value would still get passed to the user so someone might actually notice something is really wrong.  With abs(S), we'd never know.

Tim Mitchell <tmitchell>
Thu 14 Feb 2019 06:23:33 PM UTC, comment #3: 

Attached is a sample changeset that works for me.

(file #46264)

Rik <rik5>
Group administrator
Thu 14 Feb 2019 05:39:30 PM UTC, comment #2: 

Matlab doesn't display the negative sign for zero, so you need to infer it from something like taking 1/x.

Here is a sample session from R2018b showing that Matlab is returning a -0.


n = 10;
A = randn(n);
A = triu(A);

A(1) = -0;
[U,s,V] = svd(A);
1/s(end)
ans = -Inf


Also, for some reason this is dependent on the random values in A.  A more interesting test is


n = 10;
for i = 1:20
  A = randn(n);
  A = triu(A);
  A(1) = -0;
  [U,S,V] = svd(A);
  tst(i) = (1/S(end) == -Inf);
end
fprintf ("-0 results: %d / 20\n", sum (tst));


This shows about half the time the results are -0 versus +0.

This occurs on both Matlab and Octave.

I'm pretty sure this will be found to be an upstream issue with LAPACK, but in the meantime, we can probably do something about it.

The matrix S returned from the [U,S,V] is a diagonal matrix which in Octave means it is stored in a space-efficient manner.  If we were to call abs(S) before returning it would only require running abs() on the N elements of the diagonal which would be exceedlingly quick.  Alternatively, we could fix things up in the m-file with


S(S == 0) = 0;


which would replace any negative zeros with positive zeros.

@Tim: Is that an acceptable workaround?  Are all other singular values guaranteed to be positive so that abs or the indexing trick would be enough?

Rik <rik5>
Group administrator
Thu 14 Feb 2019 04:04:56 PM UTC, comment #1: 

It appears to be fixed since 5.0.0.

Marco Caliari <caliari>
Group Member
Thu 14 Feb 2019 08:33:10 AM UTC, original submission:  

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).


Tim Mitchell <tmitchell>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #46264:  bug55710.cset added by rik5 (1KiB - application/octet-stream)
file #46251:  svdzero.m added by tmitchell (2KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by tmitchell (Submitted the item)
  • -email is unavailable- added by tmitchell
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-15 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2019-02-14 rik5 Attached File- Added bug55710.cset, #46264
        StatusConfirmed Patch Submitted
        Release4.0.0 dev
    2019-02-14 rik5 StatusNone Confirmed
    2019-02-14 tmitchell Attached File- Added svdzero.m, #46251
        Carbon-Copy- Added tmitchell

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code