7
$\begingroup$

Update 2  Per a request from Daniel Lichtblau at Wolfram Research, a minimal example has been uploaded as a pure-ASCII no-format "m"-file named "oneMatrixThatFailsSVD_2017.m".

The file is rather large (about 7.7 MBytes) because it encodes a convergence-failing 320x320 complex matrix as an integer byte-array of dimension 320x320x16. This bit-perfect encoding "trick" is necessary because SVD convergence-failure seemingly is exquisitely sensitive to the least-significant bits of IEEE complex numbers (which are of course precisely the bits to which no well-conditioned SVD algorithm should be sensitive).


Update  On Mar 11, 2016, near-identical symptoms (namely, SVD convergence-failure) were reported on GitHub's "Anaconda-Issues" forum as "BUG: LinAlgError raised by scipy.linalg.svd for valid matrix #695".

By Apr 26, 2016, Anaconda users reported that:

Great news -- the Intel folks could replicate the issue with the same array 
that Anaconda currently fails on, and they report that it is fixed 
by the current pre-release version of MKL 11.3.3

Here "MKL" is an acronym for (afaict) the "Math Kernel Library".

However, as of March 9, 2017, Intel users and engineers were still wresting with this problem, as reported in the thread "Bug in GESDD (but not GESVD)". It is not apparent that the source of the SVD failures is understood, much less that a software fix has been implemented.

Conclusion  Multiple lines of evidence seeming indicate that

  1. an SVD bug exists,
  2. that is reproducibly elicited by unremarkable matrices,
  3. that has been affecting multiple users,
  4. on multiple systems (including Mathematica),
  5. for multiple years,
  6. whose algorithmic origin(s) are unknown,
  7. and for which no software fix(es) are presently available.

The bug in a nutshell:  for "MachinePrecision" complex matrices, Mathematica's SingularValueDecomposition[] sporadically fails to converge.

The associated Mathematica-generated error message is:

SingularValueDecomposition::cflsvd: Machine-precision algorithm
  failed to converge. Arbitrary-precision algorithm is called, 
  which is slower but more accurate.

This is a followup to a long-standing Mathematica bug report (specifically Wolfram Research bug report [TS 28968], submitted way back in 2005).

This tarball provides (in a folder named "SVDfailures_2017") 25 examples of matrices whose convergence fails under Mathematica 10.2.0 for Mac OS X x86. The same tarball provides (in a folder named "SVDfailures_2005") matrices that fail under various versions of Mathematica dating back to 2005.

A principal difference between the 2005 failures and the 2017 failures is that (at least some of) the matrices that failed outright back in 2005, now generate the (undocumented?) convergence-failure message "SingularValueDecomposition::cflsvd"

To anticipate some questions:

  • The arbitrary-precision evaluation does yield a correct decomposition, at the expense of a runtime that is 500-1000X longer.
  • The matrices that fail of convergence are (seemingly) unremarkable in respect to numerical condition and rank.
  • There is no reason (known to me at least) why SVDs of $\mathcal{O}(1)$-entry "MachinePrecision" matrices should ever "fail to converge", and there is no linear algebra software other than Mathematica's (known to me) that exhibits a similar convergence failure.

My questions are:

  • Does "SingularValueDecomposition::cflsvd" convergence-failure occur more generally, i.e., on systems other than Mathematica 10.2.0 for Mac OS X?
  • What's the best way to report this bug (if it is a reproducible bug)?
$\endgroup$
15
  • 2
    $\begingroup$ Why single-precision? Mathematica uses double-precision. $\endgroup$
    – Carl Woll
    Apr 22, 2017 at 0:19
  • 1
    $\begingroup$ My laptop does not have whatever software is needed to contend with a tarball of .mx files. It can extract them but does not tell me the location and does not understand to open them using Mathematica. The effort to address this got my machine into some deadlocked state requiring a reboot. I am simply out of patience. If you want an example inspected, make a link to an individual .m file containing one sample matrix. It is around 3 Mb per matrix, I believe, so this should not require a tarball. $\endgroup$ Apr 22, 2017 at 17:54
  • 2
    $\begingroup$ Link to Wolfram Community thread on same topic. $\endgroup$ Apr 22, 2017 at 19:57
  • 2
    $\begingroup$ You can get full 16 or 17 digit output using Put. I gather this is not enough, and you require actual bit patterns preserved? If so, this raises the question of how such matrices were generated: it might be simpler just to provide code to generate one. $\endgroup$ Apr 22, 2017 at 22:22
  • 1
    $\begingroup$ I can open the .m in a notebook but the fine control required to read the unpacking instructions t the bottom is not working. That is to say, I can scroll but I cannot use the line-wise up and down arrows to get where I need to be (that is, I cannot see that "Comment at the end..."). If instead I use Get[] on the .m file, it runs all the code, which I do not want, and I have no idea what the matrix is called. If I try to define it by typing "mat=" at the top of the notebook, it takes quite a long time to evaluate. Might work, maybe. I'll try at work tomorrow, opening in a text editor first. $\endgroup$ Apr 23, 2017 at 16:24

1 Answer 1

6
$\begingroup$

I do not have an answer but I want to avoid the "avoid extended discussion in comments" and anyway this is too long for comments.

First: I do get that same failure to converge message when I run the code in the linked .m file. I do not know what particular route it takes in LAPACK so I'm not sure how to track it. I might be able to get a bit of insight when I try tracking the route of the bignum version, if it happens that our bignum code accurately mirrors the LAPACK version.

Next: "...SVD convergence-failure seemingly is exquisitely sensitive to the least-significant bits of IEEE complex numbers (which are of course precisely the bits to which no well-conditioned SVD algorithm should be sensitive)." This is true but it also provides a hint as to what might be the general nature of the problem. Also I did some perambulation through LAPACK sources.

http://www.netlib.org/lapack/explore-html/index.html

http://www.netlib.org/lapack/explore-html/d0/da6/group__complex16_o_t_h_e_rcomputational_ga42f492d0f5e62a073a80f9ae57a5ee62.html#ga42f492d0f5e62a073a80f9ae57a5ee62

Comments in those files suggest that there may be (uncommon) failure-to-converge states. What happens I suspect has to do with something as basic as trichotomy: the assumption that for real numbers (a,b), either a<b or a=b or a>b. The fact that manifestations depend on machine epsilonic differences in input values, and might even be processor dependent (per discussion in one of the links), are what incline me in this direction.

One way to get a failure is to have an actual mismatch in the code, wherein one place uses less-than and another uses less-equal. This can led to a loop that in effect starts to do nothing at some point and continues this for all successive steps but fails to recognize it has achieved convergence. The people who write and maintain LAPACK know what they are doing and most likely did not make this mistake. That said, given the thicket of code paths it is always possible that such a mismatch exists between some particular pair of routines that use different conventions for checking convergence e.g. IF( thetamin .LT. thresh ) vs. IF( thetamin .LTE. thresh ).

Another way to run afoul of this, somewhat more subtle, is to have a pair of values that explicitly violates trichotomy. Of course that's always impossible...except when it's not. Which is to say, I myself ran afoul of this a couple of years ago. The way it happens, in machine arithmetic, is when some code twice computes the same numeric value, possible but not necessarily in two different ways, and then compares them. They should be equal. Depending on vagaries such as optimization level, vector 8 vs 16 byte alignment, and maybe other details, one value might or might not remain in a register that, depending on architecture, might or might not be larger than a machine double (typical is 80 bits vs 64 for machine doubles). If a computed value has not been stored in a 64 bit location but instead kept in its register, a comparison might well determine that it is strictly greater, or less, than the same value computed earlier and stored in memory.

Such a comparison error can then lead to a loop wherein one part thinks no more work is needed to get convergence, and another thinks convergence has not been attained. Is this the cause of the problems in this particular case? Obviously I don't know. All I can say is this general type of problem is consistent with all the data I have seen that describes circumstances under which is is seen.

--- edit ---

I got the .m to behave for me and it seems to generate the correct complex matrix. Steps:

(1) Open the .m file in a text editor.

(2) Precede the 320x320x16 integer array with the assignment: anIntegerArrayThatFails = ...

(3) Comment out everything after that definition.

In Mathematica kernel we now proceed as follows.

(4) Load the .m file. I just use Get[] for this.

(5) Recreate the matrix:

cmat = Partition[ ImportString[ExportString[anIntegerArrayThatFails, "Byte"], "Complex128"], Length[anIntegerArrayThatFails]];

(6) Sanity checks:

In[4]:=Dimensions[cmat]

Out[4]= {320, 320}

In[5]:=InputForm[cmat[[10,11]]]

Out[5]//InputForm= 0.0005382381370357449 - 0.00007306388005845672*I

(6) Now do the SVD:

In[6]:= svd = SingularValueDecomposition[cmat];

For better or worse it works fine, on the Linux machine I ran it on in the current Mathematica kernel. I show the singular values below.

In[7]:= Diagonal[svd[[2]]]//InputForm

Out[7]//InputForm= {2.904828261852695, 2.2951344541550434, 1.0918227589769098, 1.004862161464147, 0.8610225469711239, 0.6681124564398183, 0.6349682092830861, 0.6103146300444966, 0.5680816716691716, 0.5426987022981051, 0.49489169787986537, 0.46129735194393734, 0.43333194698189315, 0.41842507370177495, 0.41416834742927083, 0.3945741797105204, 0.38883686437150744, 0.3804953135414925, 0.36091214263408594, 0.3360101799616245, 0.3236909086672861, 0.3143821092134272, 0.3049002562655973, 0.30010582742634107, 0.29123291234527315, 0.28027130047535653, 0.26477078882783617, 0.236967943806251, 0.22557102384255504, 0.2102944648644782, 0.19238550015013525, 0.1826389902666694, 0.17336157778893144, 0.16445221904752633, 0.15731555318493404, 0.15470127447783716, 0.1456334110900463, 0.14280464903626908, 0.1381578565960464, 0.1344236796378875, 0.1310993498361891, 0.12666815647613922, 0.12605974491804625, 0.11969335869485165, 0.11895935616324389, 0.11707128220546593, 0.11420264975737747, 0.1103268165070405, 0.10738921091580042, 0.10671367357149583, 0.10367467681566253, 0.10195791921075832, 0.09974277876197389, 0.0983402105969247, 0.09589818665400816, 0.09394400698833252, 0.09207289847672227, 0.09012811915969542, 0.08829470819544377, 0.08619036258636893, 0.08332547501008657, 0.08270298910596832, 0.08042448719860448, 0.0789055591619517, 0.0779663868877628, 0.07622547020175933, 0.07526331360903826, 0.0739398832516949, 0.07284387542304292, 0.07019099016301632, 0.0678054341289208, 0.06553471311074177, 0.06454946492413614, 0.06289139925798176, 0.061011396806530435, 0.0591993281158084, 0.05794770963746927, 0.056204757179840775, 0.05515527237268259, 0.05372875315183203, 0.051206427560817766, 0.049203326736377014, 0.0479307547432323, 0.04712585208152816, 0.04593727104009871, 0.0446884313420753, 0.04375696699226062, 0.04201504537678756, 0.041465841221084404, 0.04033114255574, 0.0377828945967318, 0.0370894102987785, 0.035677653876427205, 0.033425182604297714, 0.03279900861658223, 0.031570514826654045, 0.029763840078056197, 0.029440915559297048, 0.028961728414516608, 0.028066820494299965, 0.02710620440869419, 0.026596330242566763, 0.02525814101350337, 0.02445795199731368, 0.023650534716997223, 0.02280962983108782, 0.02204232918048599, 0.021391422191204705, 0.02097280337529124, 0.02010277066977184, 0.019676786987286563, 0.018921687651602313, 0.01838391130180509, 0.017742161837301127, 0.017130206522529806, 0.016891101053555808, 0.01624292899184672, 0.015765427185192273, 0.01494655583960871, 0.01441425711006238, 0.014083731482745056, 0.013575267318266379, 0.013008938102496814, 0.012458409857924747, 0.012089734437334815, 0.011583960591848247, 0.010949282041555706, 0.010578492368751616, 0.010418547015113123, 0.009758914398614923, 0.009266255375021137, 0.009142257351872423, 0.008939517550163324, 0.008403994805638507, 0.008210067339576506, 0.0076438295086334405, 0.007339445540626259, 0.006955252203215465, 0.006705436226624429, 0.006473784884741475, 0.00636003007790285, 0.005988070077346958, 0.005805681341589245, 0.005541562510167515, 0.005176422394760754, 0.004830746770307958, 0.004736367580889274, 0.004540494242104383, 0.0043766362115849325, 0.004141239429134337, 0.003928204950297248, 0.0035823149039006844, 0.0033354897920324735, 0.0032208796636260316, 0.003043407612025262, 0.0027927309774650045, 0.0024658310023512662, 0.00240272093832309, 0.0023367556988262935, 0.002166216049750861, 0.002031068082098831, 0.00177494276158353, 0.0016764955201414582, 0.0015858402609002864, 0.0015225078635908546, 0.0014155710312122587, 0.0012917305878034438, 0.0011850471053676474, 0.0010126251411202998, 0.000957750874990038, 0.0008348499682612247, 0.0007229707037671126, 0.0006379547294593806, 0.0005563952281300264, 0.0005232025634655019, 0.0004391439440147761, 0.0003567129094759658, 0.0003150158742233563, 0.00022492723289999745, 0.00017918736178524283, 5.973543307436641*^-13, 5.95596864761285*^-13, 5.671572294842507*^-13, 5.649692543251631*^-13, 5.628096655868066*^-13, 5.506937816485377*^-13, 5.347542265503219*^-13, 5.325265205927937*^-13, 5.219366564787722*^-13, 5.171224336783615*^-13, 5.11070998274312*^-13, 5.059898018754027*^-13, 5.006253451154594*^-13, 4.995315078572324*^-13, 4.88742732922493*^-13, 4.825012104269872*^-13, 4.778887955132331*^-13, 4.73009882323133*^-13, 4.718766868644895*^-13, 4.617691808782169*^-13, 4.615029942187599*^-13, 4.5245925039243575*^-13, 4.475253976497849*^-13, 4.419811229481811*^-13, 4.3620504858309034*^-13, 4.2465671641870307*^-13, 4.2072871069100844*^-13, 4.1550555993891344*^-13, 4.115982063335738*^-13, 4.109070454714896*^-13, 4.024748504059615*^-13, 3.9695511373757535*^-13, 3.9530267719911067*^-13, 3.8911586406379883*^-13, 3.838602448037583*^-13, 3.818218373892585*^-13, 3.7708310291417817*^-13, 3.726424587436366*^-13, 3.684649253959855*^-13, 3.6823125186240376*^-13, 3.649764271942159*^-13, 3.5880323452485686*^-13, 3.538314340662835*^-13, 3.4653460661908166*^-13, 3.4591447206878043*^-13, 3.3800533290956435*^-13, 3.3306779823303953*^-13, 3.325343531319368*^-13, 3.292292390761577*^-13, 3.209506655124878*^-13, 3.1841411852563575*^-13, 3.183627098191257*^-13, 3.110835178398297*^-13, 3.074843763611202*^-13, 3.0490629799765705*^-13, 2.9670441944166715*^-13, 2.938475400835273*^-13, 2.9241263854851885*^-13, 2.921901029037186*^-13, 2.8544041350432876*^-13, 2.8033571465846197*^-13, 2.759964271001196*^-13, 2.6741190683714956*^-13, 2.638371662454207*^-13, 2.635816235341088*^-13, 2.5641521269828595*^-13, 2.5603337639199134*^-13, 2.50921206845509*^-13, 2.460450279084003*^-13, 2.459460622624308*^-13, 2.4173783692019986*^-13, 2.416945897687797*^-13, 2.317394892821592*^-13, 2.3070502038520817*^-13, 2.301170231662558*^-13, 2.2474428335551227*^-13, 2.2423682430888825*^-13, 2.1755053150774122*^-13, 2.1701335781705964*^-13, 2.1523254133574016*^-13, 2.089497409229387*^-13, 2.0643876649039162*^-13, 2.0200660587468303*^-13, 2.013312229992671*^-13, 1.9484598170000193*^-13, 1.9184158070098366*^-13, 1.848612906289478*^-13, 1.8353762534058395*^-13, 1.7833001596855028*^-13, 1.7716664750443545*^-13, 1.703606299899061*^-13, 1.7031937290887216*^-13, 1.6714467301096205*^-13, 1.624415968597149*^-13, 1.5980948827774032*^-13, 1.5541902394417834*^-13, 1.4671559713854797*^-13, 1.4653693567032798*^-13, 1.4289352712330218*^-13, 1.4197074073693232*^-13, 1.3619954879716156*^-13, 1.3491441452012306*^-13, 1.2787648047217908*^-13, 1.2352865275611754*^-13, 1.2270893072583402*^-13, 1.1807525657955718*^-13, 1.1740183116859696*^-13, 1.1159136234767948*^-13, 1.0769118497211401*^-13, 1.0318467771548448*^-13, 1.0317780181626016*^-13, 9.255175839736377*^-14, 9.254498809938708*^-14, 9.012992881660789*^-14, 8.645720815602*^-14, 8.337747383908351*^-14, 7.911740911170735*^-14, 7.86439618961678*^-14, 7.410603592399277*^-14, 7.187922716908003*^-14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

The values really should not matter; this should just work. So the problem does seem to be either processor or OS dependent.

--- end edit ---

--- edit #2 ---

Here is some further info:

The failure also appeared on my Windows laptop.

At the point it leaves Mathematica code it goes into Lapack's zgesdd(). It takes path 6 there and goes to among other routines dbdsdc(), which can fail in its divide-and-conquer. The main loop has FORTRAN code shown below.

    IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
  417 *
  418 *           Subproblem found. First determine its size and then
  419 *           apply divide and conquer on it.
  420 *
  421             IF( i.LT.nm1 ) THEN
  422 *
  423 *              A subproblem with E(I) small for I < NM1.
  424 *
  425                nsize = i - start + 1
  426             ELSE IF( abs( e( i ) ).GE.eps ) THEN
  427 *
  428 *              A subproblem with E(NM1) not too small but I = NM1.
  429 *
  430                nsize = n - start + 1
  431             ELSE
...

For the reason I showed earlier it is possible that line 426 evals as true, that is, we might have abs( e( i ) ).GE.eps and at the same time abs( e( i ) ).LT.eps if they are "equal" but in one case the value abs(e(i)) was in a large register but not in the other case. Could this cause a convergence failure? Probably not (but what do I know?) Further call tree parts include dlasd1, dlasd2, dlasd3. There are various failure modes in some of these. At this point it would take someone much more familiar with the Lapack code base, and with the appropriate debugging capabilities, to get any further.

--- end edit #2 ---

$\endgroup$
11
  • 1
    $\begingroup$ (1) I'll poke around some more but I can say two things now. One is that the issue is almost certainly not related to double-vs-single precision comparison tests. All code in question is likely to be uniform in which type it uses and the tests are always based on the inputs they receive (the inequality itself does not specify data types so it is up to the language as to how to handle comparisons across types)... $\endgroup$ Apr 23, 2017 at 22:09
  • 2
    $\begingroup$ (2)...The other thing is that this is not a problem I can fix. It happens in compiled library code and unless someone at WR knows how to build it from sources, and can make the problem appear in that way, there is no chance for us to actually find it. Even if we found it, we cannot fix ourselves; it would need to go to the LAPACK team. I'm sure they would be more than happy to address it. But keep in mind the steps involved in our even having a chance to locate it are not at all favorable. $\endgroup$ Apr 23, 2017 at 22:12
  • 1
    $\begingroup$ I'll file this as a suggestion report. $\endgroup$ Apr 24, 2017 at 18:57
  • 2
    $\begingroup$ (2) As for the revenue/open-source/who-owes-whom-what, we may have very different views. Mine, which I have taken into several NSF review panels, is that taxpayer funded work should be available to the entities who pay those taxes, and that includes corporations paying taxes in the US. For what it's worth, the NSF seems not at all averse to this view. $\endgroup$ Apr 25, 2017 at 16:41
  • 3
    $\begingroup$ (3) I'm not dead set against putting time into a problem like this. It is, after all, a fairly serious one and it hits many people in its badness. That said, I've put considerable hours into it already, some my employer's and some my own, and gotten about as far as I reasonably can. When it comes to dealing with library code, my shortcomings are best summarized by Dirty Harry: "A man's got to know his limitations". $\endgroup$ Apr 25, 2017 at 16:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.