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

[Tim Mitchell]

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

@count: for me, your two demos show GESVD being much more accurate than GESDD.
 Perhaps someone else here can test too.

In your first demo, you accidentally reversed s1 and s2.  On R2017b, for N=26,
it actually shows that GESVD returns the B matrix (exactly, not surprisingly)
while GESDD incurs large errors in the singular values (as in 1e14 relative
errors bad).  Here is a corrected version, which outputs the singular values
and then the relative errors: 


function demo1(N)
  s = logspace(1,-30, N)';
  B = diag(s);  % diagonal matrix
  
  %svd_driver('gesvd');          
  [s1] = svd(B);         % GESVD always
  %svd_driver('gesdd');           
  [u, s2, v] = svd(B);   % GESDD for N >= 26 on MATLAB
  s2 = diag(s2);
  
  % Compare the singular values:
  % truth, GESVD, GESDD (GESVD if N <= 25)
  [s s1 s2]     
  
  % The relative errors of the singular values of
  % GESVD (left) and GESDD (left)
  % For N <= 25, both will be GESVD on MATLAB
  [(s1-s)./s (s2-s)./s ]
end


On R2017b, I can't replicate the bad results you report for N=11.  For N=26,
GESVD computes each singular value to at least 11 digits while GESDD has
relative errors up to 10e23.  Here is the updated code, with same output
updates I made for demo1:


function demo2(N)
  % Make a orthonomal matrix.
  % Use qr() to avoid the doubtful svd() used in orth().
  [A0, ~]   = qr(rand(N));
  s         = logspace(0,-40,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 s1 s2]
  
  % The relative errors of the singular values of
  % GESVD (left) and GESDD (left)
  % For N <= 25, both will be GESVD on MATLAB
  [(s1-s)./s (s2-s)./s ]
end


In all the cases and tests I've seen so far, GESVD is much better than GESDD
in terms of relative accuracy of the singular values, even though it too is
not perfect.  See Figs. 29 and 30 of the recent SIAM Review article on the SVD
by Dongarra et al <https://doi.org/10.1137/17M1117732>, which precisely
demonstrate how both can be very bad but GESDD is generally much worse
accuracy-wise on poorly scaled matrices.  Meanwhile, the figures also show how
the Jacobi and QRP alternatives do handle such poorly scaled problems to high
relative accuracy, but from Figs. 27-28, are generally slower than GESVD.  By
the way, LAPACK already has at least two Jacobi-based SVD routines: dgesvj and
dgejsv (both are one-sided Jacobi, the latter with preconditioning).
 

    _______________________________________________________

Reply to this item at:

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

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