xGESDD: vastly different singular values when vectors are also requested

[t-mitchell]

Last night I found what appears to be a bad bug when computing a full SVD using xGESDD. On the particular matrix, companion_demo(n) from EigTool with n>=26, the computed singular values are totally different than those when also requesting the singular vectors. When only the singular values are requested, both xGESDD and xGESVD seem to work fine, regardless of n. But when n>=26 and singular vectors are requested, xGESDD fails.

I actually ran into this problem using MATLAB and wrote a demo to replicate the issue. Thanks to Martin Köhler for quickly translating my demo to Fortran. Both demos, with documentation, are provided. To run the MATLAB demo, you will also need to install EigTool to generate the problematic companion matrix. Just drop the .txt extensions on the Fortran and MATLAB demos (added so I could upload them here).

Many thanks,
Tim

README.txt

svdbug.f90.txt

svdbug.m.txt

[langou]

@oamarques wrote: Confirmed for N=26 (see below) with gfortran and ifort. Note that only S(1) is the same in both cases. It looks like some kind of splitting/scaling is not working correctly...
Osni

DGESDD -- Singular Values Only -- INFO = 0
S( 1) = 0.6089013695831730E+27
S( 2) = 0.6532088430598958E+10
S( 3) = 0.1000000000000001E+01
S( 4) = 0.1000000000000000E+01
S( 5) = 0.1000000000000000E+01
S( 6) = 0.1000000000000000E+01
S( 7) = 0.1000000000000000E+01
S( 8) = 0.1000000000000000E+01
S( 9) = 0.1000000000000000E+01
S( 10) = 0.1000000000000000E+01
S( 11) = 0.1000000000000000E+01
S( 12) = 0.9999999999999998E+00
S( 13) = 0.9999999999999998E+00
S( 14) = 0.9999999999999998E+00
S( 15) = 0.9999999999999998E+00
S( 16) = 0.9999999999999998E+00
S( 17) = 0.9999999999999998E+00
S( 18) = 0.9999999999999998E+00
S( 19) = 0.9999999999999998E+00
S( 20) = 0.9999999999999998E+00
S( 21) = 0.9999999999999996E+00
S( 22) = 0.9999999999999996E+00
S( 23) = 0.9999999999999993E+00
S( 24) = 0.9999999999999993E+00
S( 25) = 0.6623400625760755E+00
S( 26) = 0.1530872408329715E-09
MIN(S) = 0.1530872408329715E-09

DGESDD -- Singular Values and Vectors -- INFO = 0
S( 1) = 0.6089013695831730E+27
S( 2) = 0.6079651760069946E+11
S( 3) = 0.6079651760069946E+11
S( 4) = 0.6079651760069946E+11
S( 5) = 0.6079651760069946E+11
S( 6) = 0.6079651760069946E+11
S( 7) = 0.6079651760069946E+11
S( 8) = 0.6079651760069946E+11
S( 9) = 0.6079651760069946E+11
S( 10) = 0.6079651760069946E+11
S( 11) = 0.6079651760069946E+11
S( 12) = 0.6079651760069946E+11
S( 13) = 0.6079651760069946E+11
S( 14) = 0.6079651760069946E+11
S( 15) = 0.6079651760069946E+11
S( 16) = 0.6079651760069946E+11
S( 17) = 0.6079651760069946E+11
S( 18) = 0.6079651760069946E+11
S( 19) = 0.6079651760069946E+11
S( 20) = 0.6079651760069946E+11
S( 21) = 0.6079651760069946E+11
S( 22) = 0.6079651760069946E+11
S( 23) = 0.6079651760069946E+11
S( 24) = 0.6079651760069946E+11
S( 25) = 0.6079651760069946E+11
S( 26) = 0.2336599984275705E+10
MIN(S) = 0.2336599984275705E+10

[victorliu]

ILAENV splitting size threshold is set at 25 by default. It looks like for N <= 25 and no vectors case, the traditional algorithm is being called (DLASDQ), while for N >= 26 the divide and conquer algorithm is used (DBDSDC). Given that the smallest eigenvalues returned are about machine epsilon times the largest, is this result truly that unexpected? D&C is only known to possess high absolute accuracy in eigenvalues, not high relative accuracy.

…

On Jan 23, 2019, at 6:47 PM, langou ***@***.***> wrote: @oamarques wrote: Confirmed for N=26 (see below) with gfortran and ifort. Note that only S(1) is the same in both cases. It looks like some kind of splitting/scaling is not working correctly... Osni DGESDD -- Singular Values Only -- INFO = 0 S( 1) = 0.6089013695831730E+27 S( 2) = 0.6532088430598958E+10 S( 3) = 0.1000000000000001E+01 S( 4) = 0.1000000000000000E+01 S( 5) = 0.1000000000000000E+01 S( 6) = 0.1000000000000000E+01 S( 7) = 0.1000000000000000E+01 S( 8) = 0.1000000000000000E+01 S( 9) = 0.1000000000000000E+01 S( 10) = 0.1000000000000000E+01 S( 11) = 0.1000000000000000E+01 S( 12) = 0.9999999999999998E+00 S( 13) = 0.9999999999999998E+00 S( 14) = 0.9999999999999998E+00 S( 15) = 0.9999999999999998E+00 S( 16) = 0.9999999999999998E+00 S( 17) = 0.9999999999999998E+00 S( 18) = 0.9999999999999998E+00 S( 19) = 0.9999999999999998E+00 S( 20) = 0.9999999999999998E+00 S( 21) = 0.9999999999999996E+00 S( 22) = 0.9999999999999996E+00 S( 23) = 0.9999999999999993E+00 S( 24) = 0.9999999999999993E+00 S( 25) = 0.6623400625760755E+00 S( 26) = 0.1530872408329715E-09 MIN(S) = 0.1530872408329715E-09 DGESDD -- Singular Values and Vectors -- INFO = 0 S( 1) = 0.6089013695831730E+27 S( 2) = 0.6079651760069946E+11 S( 3) = 0.6079651760069946E+11 S( 4) = 0.6079651760069946E+11 S( 5) = 0.6079651760069946E+11 S( 6) = 0.6079651760069946E+11 S( 7) = 0.6079651760069946E+11 S( 8) = 0.6079651760069946E+11 S( 9) = 0.6079651760069946E+11 S( 10) = 0.6079651760069946E+11 S( 11) = 0.6079651760069946E+11 S( 12) = 0.6079651760069946E+11 S( 13) = 0.6079651760069946E+11 S( 14) = 0.6079651760069946E+11 S( 15) = 0.6079651760069946E+11 S( 16) = 0.6079651760069946E+11 S( 17) = 0.6079651760069946E+11 S( 18) = 0.6079651760069946E+11 S( 19) = 0.6079651760069946E+11 S( 20) = 0.6079651760069946E+11 S( 21) = 0.6079651760069946E+11 S( 22) = 0.6079651760069946E+11 S( 23) = 0.6079651760069946E+11 S( 24) = 0.6079651760069946E+11 S( 25) = 0.6079651760069946E+11 S( 26) = 0.2336599984275705E+10 MIN(S) = 0.2336599984275705E+10 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[t-mitchell]

Can you please clarify which algorithm is used when N >= 26 but no vectors are requested? In that case, the singular values seem to be computed correctly. This would help shed light on whether it is really a bug or a deficiency of the D&C method itself.

[t-mitchell]

Martin Köhler and I just answered my own question directly above. When N >= 26 but no vectors are requested, xGESDD does not use D&C. So this behavior seems to be either a bug in D&C or just its expected behavior (albeit bad).

[t-mitchell]

Martin and I changed the LAPACK tuning parameter so xGESDD would switch to D&C for N >= 10 and the singular vectors requested. What we saw is that the singular values computed by xGESDD actually stop agreeing with those computed by xGESVD for N >= 18, with the smallest singular value not even agreeing to a single digit.

[t-mitchell]

Another update from Martin and I. We also added gsvd (MATLAB) and DGGSVD to the demos, to compare to the generalized singular values of A,I, which should agree. In this case, we see that the computed values stop agreeing with the regular SVD routines at N >= 18. I guess that is when DGGSVD automatically switches to D&C? Attached are our updated demos.

svdbug_v2.f90.txt

svdbug_V2.m.txt