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 ---
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$