6
$\begingroup$

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:

LinearAlgebra`MatrixConditionNumber[mat]
(*6.85852*)

LinearAlgebra`MatrixConditionNumber[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.

enter image description here


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]}]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ I cannot explain this behavior but I will note that using SparseArray @ Chop @ mat eliminates the message and gives the same MatrixConditionNumber as mat. $\endgroup$ – Mr.Wizard Mar 28 '17 at 16:26
  • $\begingroup$ It certainly is badly conditioned: i.stack.imgur.com/H1Gv6.png -- I know hardly anything about sparse solvers. $\endgroup$ – Michael E2 Mar 28 '17 at 17:04
  • $\begingroup$ @MichaelE2 According to my test, the LinearSolve[mat][pts] gives the correct results. Please see Update:) $\endgroup$ – xyz Mar 29 '17 at 4:17
  • 1
    $\begingroup$ @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. $\endgroup$ – Michael E2 Mar 29 '17 at 10:34
  • 4
    $\begingroup$ A bug report was filed on this. $\endgroup$ – Daniel Lichtblau Mar 29 '17 at 16:25
4
$\begingroup$

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)"] &

Mathematica graphics

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)"] &

Mathematica graphics

$\endgroup$
  • $\begingroup$ 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.:) $\endgroup$ – xyz Mar 30 '17 at 13:58
  • $\begingroup$ In addition, my workaround is very simple, like mat = Normal[matSparse]. $\endgroup$ – xyz Mar 30 '17 at 14:01
5
$\begingroup$

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:

histogram

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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