bugGNU Octave - Bugs: bug #55564, GESDD: vastly different singular...

 
 

bug #55564: GESDD: vastly different singular values when vectors are also requested

Submitter:  None
Submitted:  Thu 24 Jan 2019 10:03:53 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Confirmed Assigned to:  None
Originator Name:  Tim Mitchell Originator Email:  -email is unavailable-
Open/Closed:  * Open Release:  * dev
Operating System:  * BSD Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 24 Feb 2019 03:03:32 PM UTC, comment #22: 

Just in case anyone here is interested in help testing the svd_driver GEJSV, see #55727#comment2 for the patch.

It can generate reference answers to the test cases here.

 Eddy <count>
Fri 15 Feb 2019 06:56:15 PM UTC, comment #21: 

The point is: GESVD is also returning wrong result for the original example. i.e. you only gain extra precision of one or two digits after the switch to GESVD. see (4) below, and my uploaded file #46268.

@Tim Mitchell

1. Sorry for my misleading words. Yes, I mean I'm suggesting to use GESDD as the default driver, and keep GESVD as an alternative, and a proper documentation in svd().

  The reason is, as I said and as you understood: (i) error=1e-15 is not too bad; (ii) GESVD too slow for large scale problem.

  The speed is alread discussed elsewhere, but let me post one more sample.

a = randn(2000); tic; [u,s,v]=svd(a); toc;
octave, GESVD: 54.822 sec
octave, GESDD: 3.6 sec
matlab: 3.2 sec

  We gain extra (say 2 digits) precision occasionally at the cost of 10 times slower, every time, despite of the user need it or not (assume average users keep using the default driver). Most users will feel that Octave is slow, instead of been more accurate.

2. For the dgesvj or dgejsv drivers. I will file a new bug.

3. Thanks for linking the SIAM Review. It helps a lot. Thanks for correcting the relative error computations.

4. Right, GESDD returns wrong result there, but GESVD is also wrong for the same example. The second singular value of A:


N = 26
+nomarkup+ formula solution = 1
+nomarkup+ GESVD = 3.76175e+09
+nomarkup+ GESDD = 6.07965e+10

Clearly, this example can't be used to demonstrate that GESVD is better than GESDD. They are all wrong, well GESVD is less wrong by one digit.

Compared to the first singular value s_max = 6.08901e+26, they are all around or under machine epsilon.

In the SIAM review, fig.29 and fig.30 compared the accuracy for us. The improvement of QR (GESVD) compare to D&C (GESDD) range from zero to several digits. I think it is not worth to switch for this small and uncertain improvement.

 Eddy <count>
Fri 15 Feb 2019 04:34:36 PM UTC, comment #20: 

@count: Go ahead and file a different bug report about adding additional drivers to svd_drivers() as a feature request.  If you have software skills it would also be good to submit a patch.  This is a relatively esoteric request which may not attract anyone else to write such a patch.

Beyond that, no algorithm is perfect.  Octave is prioritizing accuracy over speed, but--unlike Matlab--we give you the choice of using whatever driver you want.  I think the documentation, quoted below, is clear on that.


 -- VAL = svd_driver ()
 -- OLD_VAL = svd_driver (NEW_VAL)
 -- svd_driver (NEW_VAL, "local")
     Query or set the underlying LAPACK driver used by 'svd'.

     Currently recognized values are "gesdd" and "gesvd".  The default
     is "gesvd".

     When called from inside a function with the "local" option, the
     variable is changed locally for the function and any subroutines it
     calls.  The original variable value is restored when exiting the
     function.

     Algorithm Notes: The LAPACK library provides two routines for
     calculating the full singular value decomposition (left and right
     singular matrices as well as singular values).  When calculating
     just the singular values the following discussion is not relevant.

     The newer 'gesdd' routine is based on a Divide-and-Conquer
     algorithm that is 5X faster than the alternative 'gesvd', which is
     based on QR factorization.  However, the new algorithm can use
     significantly more memory.  For an MxN input matrix the memory
     usage is of order O(min(M,N) ^ 2), whereas the alternative is of
     order O(max(M,N)).

     Beyond speed and memory issues, there have been instances where
     some input matrices were not accurately decomposed by 'gesdd'.  See
     currently active bug <https://savannah.gnu.org/bugs/?55564>.  Until
     these accuracy issues are resolved in a new version of the LAPACK
     library, the default driver in Octave has been set to "gesvd".




 Rik <rik5>
Group administrator
Fri 15 Feb 2019 01:52:33 PM UTC, comment #19: 

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

 Tim Mitchell <tmitchell>
Fri 15 Feb 2019 01:11:46 PM UTC, comment #18: 

...In comment #17, the typo fix is for my #15

@Tim Mitchell

For demo1(), though 'gdesdd' start to fail around 1e-15. I don't feel too bad about it. It is around machine epsilon.

For demo2(), don't change the order of the logspace(), my code is for demo that 'dgesvd' is not better than 'dgesdd' too much.

It should looks like this (note the logspace() range and the transpose of s):



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)
  sort([s s1 s2], 'descend')

  % 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


My Matlab version is 9.3.0.713579 (R2017b) under Linux, using CPU i7-6700HQ (with avx2 and fma). No idea why can't be reproduced.

 Eddy <count>
Fri 15 Feb 2019 12:50:02 PM UTC, comment #17: 

I forgot to mention that I vote for keep using the 'gesdd' driver.

There is one more SVD solver in LAPACK: 'dgejsv', which uses Jacobi SVD algorithm. It should be the most accurate one in LAPACK. If possible, I wish to add that as an addition svd_driver. Should I open another ticket/issue for this?


Fix for some typos in comment #16 ...

"demo 'gesdd' is better than 'gesvd'"
it is a demo of 'gesvd' is better than 'gesdd':

In the attachement (file #46268), The comment

% Characteristic polynomial of A*A', note svd(A)=eig(A*A')=roots(c).

should be

% Characteristic polynomial of I-A*A', note svd(A)^2=eig(A*A')=1-roots(c).

 Eddy <count>
Fri 15 Feb 2019 12:27:36 PM UTC, comment #16: 

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

Tim Mitchell <tmitchell>
Fri 15 Feb 2019 12:47:56 AM UTC, comment #15: 

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), but I don't have a good code in hand  that can handle any of above extreme case.


(file #46268)

 Eddy <count>
Fri 25 Jan 2019 07:00:08 PM UTC, comment #14: 

I re-opened bug #49940 and linked it to this report.

 Rik <rik5>
Group administrator
Fri 25 Jan 2019 06:54:27 PM UTC, comment #13: 

I made "gesvd" the default driver in Octave (https://hg.savannah.gnu.org/hgweb/octave/rev/1212568010a8).  I updated the documentation to explain the issue and with a URL reference to this bug report.

 Rik <rik5>
Group administrator
Fri 25 Jan 2019 06:30:14 PM UTC, comment #12: 

@Mike: Given that the inaccuracies of GESDD's results can be arbitrarily bad, I'd personally be wary to add any documentation encouraging users to try it, even if speed is valued more than accuracy.   If anything, until this is better resolved on the LAPACK side, I think it would be prudent to warn users that using GESDD could be trading off all accuracy for speed!

 Tim Mitchell <tmitchell>
Fri 25 Jan 2019 06:20:42 PM UTC, comment #11: 

@Mike: At the moment, no other solution than documentation.  We could look to see if there is anything within Octave, such as unnecessarily duplicating matrices, that we might tweak.  But I believe the difference in performance is really located with LAPACK.

 Rik <rik5>
Group administrator
Fri 25 Jan 2019 06:02:11 PM UTC, comment #10: 

From my understanding, bug #49940 is about the performance of the GESVD function, which we are now considering making the default again. Is there any solution for that other than documenting it and suggesting that users who care more about speed than accuracy should use GESDD?

 Mike Miller <mtmiller>
Group Member
Fri 25 Jan 2019 04:17:53 PM UTC, comment #9: 

I generally agree. I guess my point was that we switched to gesdd
to resolve bug 49940. So if we go back (and I agree that at this
moment we should), we need to re-open 49940.
Also the benchmark ratio seems got worse.

Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Fri 25 Jan 2019 04:07:33 PM UTC, comment #8: 

We can look at bug #49940 again.  But, I still believe we need to prioritize correctness over speed.  The issue with this bug is that most people will be unaware that they have incorrect results.  Even the pessimists, who verify the results using the acceptance test in the documentation, will believe that the results are good.  It is only if you calculate using GESVD that you will find the problem.  Since any new acceptance test is really going to require running GESVD, we might as well do that in the first place.

I re-read the discussion in bug #49940, and I'm pretty sure it isn't Octave's problem.  My guess is that the GESDD "Divide and Conquer" algorithm is able to spread the problem over multiple threads/cores, but arrives at the wrong answer.  And the GESVD algorithm is probably old and single-threaded with few opportunities for parallelism.

 Rik <rik5>
Group administrator
Fri 25 Jan 2019 02:20:24 AM UTC, comment #7: 

May be we should revisit bug 49940.
It looks to me that gesvd mostly runs single-thread
even with multithreaded blas.

I run a simple benchmark and watched top. With
gesdd cpuload was 400% (4 threads), with
gesvd it started ~400% but  dropped quickly to 100%



LD_PRELOAD=/usr/lib64/libopenblaso.so OMP_NUM_THREADS=4 ./run-octave

octave:1> a=randn(3000);
octave:2> tic; [u2,s2,v2] = svd (a); toc
Elapsed time is 13.4637 seconds.
octave:3> svd_driver("gesvd")
octave:4> tic; [u2,s2,v2] = svd (a); toc
Elapsed time is 247.86 seconds.


similar results with liopenblasp (pthread interface) and atlas.

I am not sure if this is octave problem (posted results in 49940
did not show such a big difference), or lapack's.

Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Fri 25 Jan 2019 01:51:52 AM UTC, comment #6: 

Sure, seems prudent to me. The default was last changed in 2016 for the 4.4.0 release.

I get identical test suite results with both drivers, and confirm switching to GESVD fixes the problem reported here.

 Mike Miller <mtmiller>
Group Member
Fri 25 Jan 2019 01:09:32 AM UTC, comment #5: 

@Mike: I think we should favor correctness over speed.  What do you think about changing the default driver in Octave back to GESVD for the 5.1 release?

 Rik <rik5>
Group administrator
Fri 25 Jan 2019 01:06:22 AM UTC, comment #4: 

@Mike: Thanks.  I renamed the function to tst_svd on my platform.  Here is a sample session which shows changing the backend SVD function to GESVD fixes the issue.


octave:6> svd_driver ('gesdd')
octave:7> tst_svd (25)
octave:8> tst_svd (26)
error: called from
    assert at line 424 column 7
    tst_svd at line 6 column 3
error: ASSERT errors for:  assert (s1,s2,-10 * eps)

  Location  |  Observed  |  Expected  |  Reason
    (2)      6532088430.599 60796517600.6995   Rel err 0.89256 exceeds tol 2.2204e-15 by 0.9
    (3)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (4)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (5)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (6)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (7)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (8)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (9)           1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (10)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (11)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (12)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (13)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (14)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (15)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (16)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (17)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (18)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (19)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (20)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (21)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (22)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (23)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (24)          1       60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (25)     0.66234006257608 60796517600.6995   Rel err 1 exceeds tol 2.2204e-15 by 1
    (26)     1.5308724083297e-10 2336599984.2757   Rel err 1 exceeds tol 2.2204e-15 by 1
octave:8> svd_driver ('gesvd')
octave:9> tst_svd (25)
octave:10> tst_svd (26)
octave:11> diary off



 Rik <rik5>
Group administrator
Fri 25 Jan 2019 12:05:55 AM UTC, comment #3: 

Copy/paste error, this should be a little clearer:


function testit (N)
  A = compan (fliplr ([1, 1 ./ cumprod(1:N)]));
  s1 = svd (A);
  [U, S, V] = svd (A);
  s2 = diag (S);
  assert (s1, s2, -10*eps);
endfunction


testit (25)  ## ok
testit (26)  ## not ok


 Mike Miller <mtmiller>
Group Member
Fri 25 Jan 2019 12:03:00 AM UTC, comment #2: 

Here is what I've distilled out of the original problem statement on https://github.com/Reference-LAPACK/lapack/issues/316:


## N ≤ 25 works ok
N = 25;
A = compan (fliplr ([1, 1 ./ cumprod(1:N)]));
s1 = svd (A);
[U, S, V] = svd (A);
s2 = diag (S);
s1
s2

## N ≥ 26 produces different singular values
N = 25;
A = compan (fliplr ([1, 1 ./ cumprod(1:N)]));
s1 = svd (A);
[U, S, V] = svd (A);
s2 = diag (S);
s1
s2


 Mike Miller <mtmiller>
Group Member
Thu 24 Jan 2019 11:50:05 PM UTC, comment #1: 

Is the original problem matrix available in a file (.mat or otherwise)?  Could that be uploaded?

 Rik <rik5>
Group administrator
Thu 24 Jan 2019 10:03:53 PM UTC, original submission:  

I recently opened an issue on LAPACK's Github page about how the GESDD routine can produce vastly different singular values if the singular vectors are also requested.  At the moment, it is still unclear if this behavior is a coding bug or just an inherent deficiency of the "divide and conquer" approach that GESDD implements.  You will need to follow the LAPACK issue for developments on that and for complete details.

As Octave 4.4 is using GESDD by default (like MATLAB and others), I just wanted to raise your attention of how bad this problem can be in practice.  For future Octave releases, you might consider only using GESVD for all svd calls, unless the user specifically enables GESDD (at their own risk). ;-)

In any case, GESDD's inaccuracy in singular values can basically be arbitrarily inaccurate.  To be fair, note that this may only happen rarely and indeed the example I found was designed to be numerically challenging (companion_demo(n) from EigTool, with n>=26).  Nevertheless, even though this matrix is nearly singular (GESVD returns 2.6e-10 as the smallest singular value), GESDD instead returns 2.3e+9 as the smallest singular value!  For this problem, only the largest singular value is computed consistently; all others don't even agree to a single digit between GESVD and GESDD. 

Lastly, note that using norm(U*S*V'-A)/norm(A) to verify the computation for this example is not a sufficient test.  The reason is that the largest singular value is 6e+25 and all the others are at least 1e+16 times smaller.  So, to machine precision, the matrix seems to be rank one.  In other words, one can randomly choose any value in [0,1e9] for the other n-1 singular values and still have norm(U*S*V'-A)/norm(A) look good numerically!  A better test is to compare the singular values to those computed by GESVD and make sure each singular value has a high number of digits in agreement.

 Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46268:  svd_inacc2.m added by count (2KiB - text/x-objcsrc - Test relative accuracy of svd)

 

Depends on the following items: None found

Items that depend on this one
bugsThe svd function takes very long time for big matrices (bug #49940, Confirmed)

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by count (Updated the item)
  • -email is unavailable- added by count
  • -email is unavailable- added by tmitchell (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)

    •  

    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 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-02-26 mtmiller Release5.0.1 dev
    2019-02-15 count Attached File- Added svd_inacc2.m, #46268
        Carbon-Copy- Added count
    2019-01-25 rik5 Dependencies- bugs #49940 is dependent
    2019-01-25 mtmiller StatusNone Confirmed

    Back to the top

    Copyright © 2024 Free Software Foundation, Inc.
    Verbatim copying and distribution of this entire article is permitted in any medium, provided this notice is preserved.
    The Levitating, Meditating, Flute-playing Gnu logo is a GNU GPL'ed image provided by the Nevrax Design Team.
    Page source code

    Powered by Savane 3.13-7a7b.
    Corresponding source code