[Octave-bug-tracker] [bug #55564] GESDD: vastly different singular values when vectors are also requested

[Tim Mitchell]

Follow-up Comment #19, bug #55564 (project octave):

@count: 

1. I don't think anyone is proposing to remove any of the different SVD
drivers.  We are just discussing to not use GESDD by default (because of its
worse accuracy than GESVD).  

2. Yes, in my #15, I already mentioned that both dgesvj and dgejsv are
available in LAPACK.  Perhaps it would be nice if Octave could also add
support for them to be chosen as other alternatives if the user wishes.

3. Regarding demo2(), sorry about changing the order of the log space and yes,
now I can replicate that.  However, again no one was saying GESVD is perfect,
just that it is much better than GESDD; again, see the SIAM Review article. 
Also, your update completely broke the relative error computations.  Below is
a corrected version.  For N=26, while both methods indeed struggle with the
small singular values, GESDD is still *much* worse (1e23 relative error versus
1e11).  Such an example would be a reason to try the Jacobi routines.


function demo2(N)
  % Make a orthonomal matrix.
  % Use qr() to avoid the doubtful svd() used in orth().
  [A0, ~]   = qr(rand(N));
  s         = logspace(-40,0,N)';
  C         = s .* A0;
  
  s1        = svd(C);   % GESVD always
  [U,S,V]   = svd(C);   % GESDD for N >= 26 on MATLAB
  s2        = diag(S); 
  
  % Compare the singular values: 
  % truth, GESVD, GESDD (GESVD if N <= 25)
  s         = sort(s,'descend');  % Need to resort s for here and below!
  [s s1 s2]
  
  % The relative errors of the singular values of
  % GESVD (left) and GESDD (right)
  % For N <= 25, both will be GESVD on MATLAB
  [(s1-s)./s (s2-s)./s ]
end


4. GESDD does not fail at around 1e-15; it just seems to return any singular
value below *norm(A)**eps as about *norm(A)**eps, regardless of its true
value.  This means such singular values can be *completely wrong* (as in *no
digits accurate*) and this can even happen to *large* singular values if
norm(A) is very large.  That is demonstrated in the original example.  

If risking such an accuracy loss for some speedup is fine for your needs,
great.  But that tradeoff doesn't seem like a good default choice for
everyone, particularly since there are lots of applications where the higher
accuracy is needed and expected.

    _______________________________________________________

Reply to this item at:

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

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