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

[count]

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

Just wait a minute...

I found the following:

* 1) 'gesvd' did *not* give the "correct" result either, for N=26.

  I worked out an explicity formula for svd(A) in this case. So we can measure
the relative error now. See attachment for code. A point is, the second
singular value should be 1.0, but we get 3.76175e+09 from 'gesvd'.

  A straightforward example that demo 'gesdd' is better than 'gesvd':


  N = 26;
  B = diag(logspace(1,-30, N));  % diagonal matrix
  svd_driver('gesdd');           % the fast svd
  [u, s2, v] = svd(B);
  s2 = diag(s2);
  svd_driver('gesvd');           % the correct svd
  [u, s1, v] = svd(B);
  s1 = diag(s1);
  [s1, s2]                       % compare the answers

ans =

   1.0000e+01   1.0000e+01
   5.7544e-01   5.7544e-01
   3.3113e-02   3.3113e-02
   1.9055e-03   1.9055e-03
   1.0965e-04   1.0965e-04
   6.3096e-06   6.3096e-06
   3.6308e-07   3.6308e-07
   2.0893e-08   2.0893e-08
   1.2023e-09   1.2023e-09
   6.9183e-11   6.9183e-11
   3.9811e-12   3.9811e-12
   2.2909e-13   2.2909e-13
   1.3183e-14   1.3183e-14
   7.5858e-16   9.9920e-16
   4.3652e-17   9.9920e-16
   2.5119e-18   9.9920e-16
   1.4454e-19   9.9920e-16
   8.3176e-21   9.9920e-16
   4.7863e-22   9.9920e-16
   2.7542e-23   9.9920e-16
   1.5849e-24   9.9920e-16
   9.1201e-26   9.9920e-16
   5.2481e-27   9.9920e-16
   3.0200e-28   9.9920e-16
   1.7378e-29   9.9920e-16
   1.0000e-30   9.9920e-16


  Here svd(B) is trivial to compute, yet, 'gesdd' is not that satisfactory.
And indeed, we see that 'gesvd' is carefully programmed.

  However, here is a frustrating example that convince you that 'gesvd' is no
much better:


  N=11;
  % Make a orthonomal matrix.
  % Use qr() to avoid the doubtful svd() used in orth().
  [A0, ~] = qr(rand(N));

  scaling = logspace(-40,0,N);
  C = scaling' .* A0;

  % two columns should be equal
  [svd(C) sort(scaling, 'descend')']

ans =

   1.0000e+00   1.0000e+00
   1.0000e-04   1.0000e-04
   1.0000e-08   1.0000e-08
   1.0000e-12   1.0000e-12
   1.0022e-16   1.0000e-16
   3.4269e-17   1.0000e-20
   1.0142e-20   1.0000e-24
   5.8744e-21   1.0000e-28
   8.6341e-25   1.0000e-32
   3.9032e-25   1.0000e-36
   4.6972e-30   1.0000e-40


  We know exactly the SVD of C by definition: U=I, S=diag(scaling), V=A0'.
Yet, 'gesvd' is disappointed, in an idealism scene.

* 2) Matlab(2017b) also returns wrong results, the numbers (both the correct
and the wrong) are the same as octave. Indicating that Matlab is also using
'gesdd'.

* 3) The error of 'gesdd' is acceptable.

  Note that the ratio between the second and the first svd(A) is about
1.0e-16, which is coincede with the machine epsilon. Similary for svd(B) and
svd(C) above. This is the precision we can *expect*, but no more. As 'gesvd'
could also produce wrong results around 1.0e-17 (in the case of svd(C) above
and svd(A) in attachement), there is no urgent need to pursuit the occasional
extra precision.

* 4) Speed is important, given that 'gesdd' is *not* terribly wrong.

So, I suggest document the difference between 'gesvd' and 'gesdd', let the
careful user pick the "correct" one.

PS: I heard that the "Jacobi" method is an even more accurate method (ref
<https://epubs.siam.org/doi/10.1137/0613074>), but I don't have a good code in
hand  that can handle any of above extreme case.


(file #46268)
    _______________________________________________________

Additional Item Attachment:

File name: svd_inacc2.m                   Size:2 KB
    <https://savannah.gnu.org/file/svd_inacc2.m?file_id=46268>



    _______________________________________________________

Reply to this item at:

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

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