# Why the sparse array is ill-conditioned?

Bug introduced in 8.0 or earlier and persisting through 11.1.0 or later

Given that I have a knot vector $$\mathbf U$$ and a set of parameter values $$\mathbf T$$:

U = {-0.018310544879457935,-0.012207034419556506,-0.0061035194597373055,0,
0.007947381024392994,0.015896985836984207,0.023849496293520733,0.0318037494975631,
0.03975915901762338,0.0477153411094262,0.05567199782586085,0.06362913932177364,
0.0715866937531541,0.0795446476499281,0.08750300546212922,0.09546173691312039,
0.10342091295646483,0.11138055246270329,0.11934063494463609,0.12730123948852812,
0.1352623862878759,0.14322404782814827,0.15118620024759044,0.15914879056899667,
0.16711176492165805,0.17507508604555258,0.18303859260454125,0.19100216524527763,
0.19896562480557684,0.20692900422432328,0.2148920912932465,0.22285488990607166,
0.23081722841868946,0.23877899625012167,0.24674016537051507,0.2547006557562886,
0.2626605712319437,0.2706197514514926,0.2785782021068993,0.28653592356892715,
0.2944929268582807,0.3024492246617148,0.3104048564483061,0.3183598587228994,
0.3263142945909529,0.3342681903482628,0.3422215816384908,0.3501747584647508,
0.35812755237990723,0.36607998995330615,0.37403190761856553,0.3819835289304483,
0.3886094255757332,0.39390964457488653,0.39788430888435367,0.4018589782161606,
0.40583364754796747,0.40966627849998927,0.41589246332954444,0.4245334531537193,
0.4357312838584109,0.4469504009218613,0.45816957544280623,0.4693888401185814,
0.4806081760209349,0.4918276023825066,0.5030471077373718,0.5142667313315674,
0.5254864642566899,0.5367063130793907,0.5479262509429471,0.5591463163437205,
0.5703663844157385,0.5815866951490725,0.5928073650583843,0.6040287621336593,
0.6152508517822124,0.6264738170541829,0.6376976658624192,0.648921694734309,
0.6601461878038705,0.6713717498203953,0.682599803804294,0.693831315188334,
0.7050674703279312,0.7163084579237454,0.7275527335316063,0.7387983309924002,
0.7500439354151093,0.7612892062041845,0.7725334514413666,0.7837761124526615,
0.7950166325928549,0.8062547781628725,0.8174905870790532,0.8287243101553067,
0.8399563055094329,0.8511870185871537,0.862416743337147,0.873645882662241,
0.884874491754384,0.8961027592682634,0.9073304502318396,0.9185573329506938,
0.9297828880613659,0.9410061871055874,0.9522046631777067,0.961696336629296,
0.9694824342007392,0.9755859446606406,0.9816894551205421,0.9877929655804435,
0.9938964805402627,1,1.007947381024393,1.0158969858369842,1.0238494962935207};

T = {0.,0.007945873121881325,0.01589626995129766,0.023848814437773633,
0.03180340449149092,0.03975902956342475,0.047715042997954475,0.055671950766899375,
0.0636289997127287,0.07158646748569283,0.07954461406104081,0.08750286140305062,
0.09546154092229622,0.10342080841401438,0.11138038953308387,0.11934045944101165,
0.12730105585981277,0.13526220316475993,0.1432238998390551,0.15118604048062978,
0.1591486604230865,0.16711167080327374,0.175074963538614,0.18303862379477004,
0.1910021904802397,0.1989656814608232,0.20692900247566773,0.21489232873647895,
0.22285494266759273,0.23081739831414336,0.23877934427433228,0.24674024616188933,
0.25470090567532366,0.2626608154316528,0.27061999258885466,0.2785784463339704,
0.28653616739787297,0.29449315697493816,0.3024494562020311,0.31040506080817526,
0.318360052334712,0.326314463025811,0.33426836841233576,0.3422217396066416,
0.3501746368964951,0.35812789889111574,0.36608012135211093,0.3740319496166919,
0.3819836518868938,0.3899349852877591,0.3939096395525468,0.3978843088843537,
0.4018589782161606,0.4058336475479675,0.4098083168797744,0.41335687107222596,
0.4245122020366331,0.43573128635229896,0.4469503631863009,0.45816955322698405,
0.469388809915134,0.4806081572136262,0.49182756093404456,0.5030470889998492,
0.5142666732782218,0.5254864317166315,0.5367062877752161,0.5479262197463246,
0.5591462453073006,0.5703664839775365,0.5815864239623786,0.5928071775073025,
0.6040284937054722,0.615250615188203,0.6264734464529619,0.6376973895213841,
0.6489221616129114,0.660145533068632,0.6713708687300681,0.6825988476624855,
0.6938296950203283,0.7050654028821882,0.7163073130812772,0.7275526578077709,
0.738798229705771,0.7500441054636585,0.7612894710758985,0.7725340420729966,
0.7837768411752049,0.7950174541097833,0.8062556024935763,0.817491277885258,
0.8287248808583252,0.839956771722337,0.8511872639476367,0.8624170200914877,
0.8736459459723168,0.8848746819229186,0.8961028473679171,0.9073307485139547,
0.9185577548136468,0.9297834955244799,0.9410074138459712,0.9522276519463112,
0.9633789237408378,0.9694824342007392,0.9755859446606406,0.9816894551205421,
0.9877929655804435,0.9938964760403449,0.9999999999999998};


Then I got a square coefficient matrix mat via corresponding computation

n = Length@T - 2;
mat = Outer[BSplineBasis[{3, U}, #2, #1] &, Most[T], Range[0, n + 3]];
{a, b} = Dimensions[mat];
mat = Take[mat, a, a] + PadRight[Take[mat, a, a - b], {a, a}];


So my confusion comes from the following operations:

LinearAlgebraMatrixConditionNumber[mat]
(*6.85852*)

LinearAlgebraMatrixConditionNumber[SparseArray@mat]
(*3.69096*10^19*)

LinearSolve[mat]

LinearSolve[SparseArray@mat]


LinearSolve::luc: Result for LinearSolve of badly conditioned matrix SparseArray[Automatic,<<3>>] may contain significant numerical errors.

### Update

Via LinearSolve[mat][pts], I can achieve a B-spline curve (defined on control points ctrlpts and knot vector U) which passes through a set of points pts.

pts = {{-15.4443,-3.0278},{-14.6621,-2.99274},{-13.88,-2.94865},{-13.0977,-2.90178},
{-12.3152,-2.85623},{-11.5324,-2.81497},{-10.7492,-2.78035},{-9.96558,-2.75406},
{-9.18168,-2.73769},{-8.39759,-2.73282},{-7.61345,-2.74044},{-6.82956,-2.76195},
{-6.04618,-2.79831},{-5.26363,-2.85047},{-4.48236,-2.91951},{-3.7028,-3.0062},
{-2.92548,-3.11148},{-2.15096,-3.23602},{-1.37985,-3.38049},{-0.612823,-3.54548},
{0.149415,-3.73153},{0.906107,-3.93913},{1.65654,-4.16843},{2.39993,-4.41973},
{3.13555,-4.69293},{3.86263,-4.98806},{4.58048,-5.30497},{5.28849,-5.64327},
{5.98603,-6.00252},{6.67259,-6.38229},{7.34777,-6.78186},{8.01107,-7.20066},
{8.66237,-7.63783},{9.30137,-8.09267},{9.92798,-8.56433},{10.5421,-9.05199},
{11.1438,-9.55477},{11.7332,-10.0718},{12.3105,-10.6024},{12.8756,-11.1456},
{13.4291,-11.7007},{13.9709,-12.2671},{14.5016,-12.8439},{15.0213,-13.4304},
{15.5305,-14.0262},{16.0294,-14.6305},{16.5183,-15.2428},{16.9977,-15.8626},
{17.468,-16.4894},{17.9292,-17.1227},{18.2309,-17.3724},{18.6154,-17.4469},
{18.9886,-17.328},{19.259,-17.0447},{19.3605,-16.6664},{19.2875,-16.3245},
{18.826,-15.3268},{18.3585,-14.325},{17.8872,-13.325},{17.4117,-12.327},
{16.9314,-11.3313},{16.4459,-10.3381},{15.9543,-9.34784},{15.4559,-8.36098},
{14.9499,-7.37804},{14.4352,-6.39961},{13.9106,-5.42643},{13.3748,-4.45933},
{12.8265,-3.49927},{12.2643,-2.54727},{11.6867,-1.60458},{11.0923,-0.672301},
{10.48,0.24846},{9.84946,1.15685},{9.19549,2.04863},{8.50955,2.9162},{7.78997,3.75618},
{7.03722,4.56638},{6.24627,5.33962},{5.41514,6.06988},{4.54293,6.75102},
{3.6291,7.37605},{2.67399,7.93719},{1.6801,8.42711},{0.651008,8.83807},
{-0.40862,9.16236},{-1.49236,9.3934},{-2.59242,9.52584},{-3.69983,9.55667},
{-4.80515,9.48539},{-5.89913,9.31365},{-6.97328,9.04539},{-8.02018,8.68584},
{-9.03385,8.24155},{-10.0094,7.71909},{-10.9431,7.12529},{-11.8323,6.46675},
{-12.6746,5.74933},{-13.4684,4.97856},{-14.212,4.15943},{-14.9039,3.29619},
{-15.5423,2.39285},{-16.1252,1.45294},{-16.6504,0.480036},{-17.1188,-0.513947},
{-17.2787,-1.09374},{-17.2344,-1.69353},{-16.9912,-2.24359},{-16.5773,-2.67996},
{-16.0409,-2.9519}};

ctrlpts = LinearSolve[mat][pts];
ctrlpts = Join[ctrlpts, ctrlpts[[1 ;; 3]]];

Graphics[{BSplineCurve[ctrlpts, SplineKnots -> U], Blue, Point[pts]}]


• I cannot explain this behavior but I will note that using SparseArray @ Chop @ mat eliminates the message and gives the same MatrixConditionNumber as mat. Mar 28 '17 at 16:26
• It certainly is badly conditioned: i.stack.imgur.com/H1Gv6.png -- I know hardly anything about sparse solvers. Mar 28 '17 at 17:04
• @MichaelE2 According to my test, the LinearSolve[mat][pts] gives the correct results. Please see Update:)
– xyz
Mar 29 '17 at 4:17
• @xyz Yes, it seemed to me that LinearSolve[mat] works fine, too, (see lsM in linked image in previous comment); but I thought it was the SparseArray one you were asking about, which is badly conditioned. Mr.Wizard's fix seems to multiply the error in lsM[b] by 1000, which is tolerable in many cases. Mar 29 '17 at 10:34
• A bug report was filed on this. Mar 29 '17 at 16:25

Here is what appears to be a slight improvement on Mr. Wizard's workaround. How much to Chop[] should probably depend on both the magnitude (norm) and precision of the matrix. Something like this:

LinearSolve[SparseArray@Chop[mat, Norm[mat] $MachineEpsilon]]  With only one test example, and given that this appears to be a bug, it's hard to test the robustness of this approach. Note that lsM = LinearSolve[mat] is accurate to machine precision. lsM = LinearSolve[mat]; b = mat.(x0 = RandomReal[1, 110]); Norm[lsM@b - #]/Norm[#] &[x0] (* 2.07362*10^-16 *)  Chopping to machine precision produces a well-conditioned matrix/linear-solve-function, whose error is only a small multiple of lsM: lsChop = LinearSolve[SparseArray@Chop[mat, Norm[mat]$MachineEpsilon]];
Table[(b = mat.(x0 = RandomReal[{-1, 1}, 110]);
Norm[lsChop@b - #]/Norm[#] &[x0])/(Norm[lsM@b - #]/Norm[#] &[x0]),
{5000}] //
Histogram[#, {0.2}, PlotRange -> All,
PlotLabel -> "Error ratio (SA@Chop[\$MachineEpsilon]/mat)"] &


The default tolerance for Chop[] has larger error (by a factor on the order of 1000):

lsChop2 = LinearSolve[SparseArray@Chop@mat];
Table[(b = mat.(x0 = RandomReal[{-1, 1}, 110]);
Norm[lsChop2@b - #]/Norm[#] &[x0])/(Norm[lsM@b - #]/Norm[#] &[x0]),
{5000}] //
Histogram[#, {100}, PlotRange -> All,
PlotLabel -> "Error ratio (SA@Chop[default]/mat)"] &


• Thanks very much for your answer and detailed analysis. In fact, I write a special function to store the non-zero values and corresponding postions, then translate them to SparseArray, rather than mat = Outer[BSplineBasis[{3, U}, #2, #1] &, Most[T], Range[0, n + 3]];. In this process, I found the above phenomenon.:)
– xyz
Mar 30 '17 at 13:58
• In addition, my workaround is very simple, like mat = Normal[matSparse].
– xyz
Mar 30 '17 at 14:01

Let me present another workaround. By default, LinearSolve[] is using a multifrontal method (Method -> "Multifrontal") on the given SparseArray[], and it seems it is having trouble with the internal thresholding, which results in the erroneous error message.

Thus, you might consider using a different method instead. In particular, the structure of the matrix involved suggests the use of Method -> "Banded":

sm = SparseArray[mat];
ls = LinearSolve[sm, Method -> "Banded"];


which gives a pretty good result in this case:

BlockRandom[SeedRandom[42];
b = sm.(x = RandomReal[1, 110]);
Norm[lsx[b] - x, ∞]/Norm[x, ∞]]
2.2219*10^-16


Here is a histogram similar to Michael's for the "Banded" strategy, with the slight replacement of using the max-norm instead of the 2-norm:

• I felt sure that someone would know or figure out how to adjust the method. Mar 30 '17 at 17:27
• It is first time for me to know that using the Method option of LinearSolve. Thanks a lot.
– xyz
Mar 31 '17 at 1:15