MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

LinearSolve[m, b] is not equivalent to LinearSolve[m][b]

  • To: mathgroup at smc.vnet.net
  • Subject: [mg84283] LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
  • From: "Andrew Moylan" <andrew.j.moylan at gmail.com>
  • Date: Mon, 17 Dec 2007 19:16:03 -0500 (EST)

At the bottom of this message I define a machine-precision matrix m and a
machine precision vector b. (I apologise for their large size.)


1. LinearSolve[m, b] is not equivalent to LinearSolve[m][b].

The documentation for LinearSolve states that "LinearSolve[m,b] is
equivalent to LinearSolve[m][b]." For m and b, however:

sol1 = LinearSolve[m, b]
>>
{-1.68129*10^9,-7.29642*10^11,-6.31285*10^20,4.99506*10^22,3.84561*10^24,-1.
3535*10^26,-7.82045*10^25,1.69903*10^27,2.25713*10^25,-2.49272*10^26,-1.0616
*10^26,9.03507*10^26,-4.77786*10^24,2.87158*10^25,4.55334*10^25,-2.43647*10^
26,-2.82578*10^26,1.27101*10^27,-1.38448*10^28,5.52681*10^30}

sol2 = LinearSolve[m][b]
>>
{-1.8893*10^-6,0.00627488,-0.000025753,-0.00767333,0.0000624517,-0.000027213
7,-0.0000464774,0.00188393,7.51337*10^-6,-0.0000585884,6.01809*10^-6,-0.0004
30135,-6.4536*10^-7,-0.0000946245,-1.29985*10^-6,0.0000701392,-4.03265*10^-7
,0.0000499027,4.97681*10^-7,9.8366*10^-6}

Needless to say, sol2 is the one that approximately solves m.x==b. Type
(m.sol1 - b) and (m.sol2 - b) to see this. sol1 is wild.


2. When using arbitrary precision arithmetic instead of machine precision
arithmetic, the wild behaviour does not occur.

Here are m and b at precision $MachinePrecision (i.e., 15.9546):
mm = SetPrecision[m, $MachinePrecision];
bb = SetPrecision[b, $MachinePrecision];

For these definitions, LinearSolve[mm,bb] and LinearSolve[mm][bb] *are*
approximately equal, and they are approximately equal to sol2 (the tame
solution).


Can anyone shed some light on this topic? My working hypothesis is that 1.
sol1 is the result of a bug (or inconsistent/bad handling of poorly
conditioned matrices); and 2. that bug isn't present in the arbitrary
precision version of LinearSolve. Then the workaround would be: don't use
LinearSolve[m, b]; instead use LinearSolve[m][b].


Andrew Moylan
 

Here are the definitions of m and b:
 
m={{0.,120.,4.,-120.,-16.,120.,36.,-120.,-64.,120.,100.,-120.,-144.,120.,196
.,-120.,-256.,120.,324.,-120.},{-120.,0.,120.,4.,-120.,-16.,120.,36.,-120.,-
64.,120.,100.,-120.,-144.,120.,196.,-120.,-256.,120.,324.},{0.,137.667770123
3155,4.,-126.52361302070159,-14.704805681953697,94.89537283269993,28.5433706
5200333,-47.90363918197348,-40.5445084603741,-6.84365962903201,45.5838278268
89575,60.48297480676686,-39.728436874545594,-104.33013932096895,21.378320417
681632,131.28632837898377,8.062099516267708,-136.987354117796,-44.1577316683
6576,120.51022451941246},{-137.6677701233155,0.,126.52361302070159,4.,-94.89
537283269993,-14.704805681953697,47.90363918197348,28.54337065200333,6.84365
962903201,-40.5445084603741,-60.48297480676686,45.583827826889575,104.330139
32096895,-39.728436874545594,-131.28632837898377,21.378320417681632,136.9873
54117796,8.062099516267708,-120.51022451941246,-44.15773166836576},{0.,190.5
9730297009568,4.,-140.79278017868444,-11.819078485136764,17.407811913964213,
14.191990544591762,115.0747432704792,-4.317884714201286,-187.4174897103089,-
15.679345960738772,161.81300927620237,34.27414014713948,-51.64259235864097,-
37.12435871470173,-85.51702774336465,16.98339066468509,177.98415019794658,19
.503939521116326,-177.43405229412292},{-190.59730297009568,0.,140.7927801786
8444,4.,-17.407811913964213,-11.819078485136764,-115.0747432704792,14.191990
544591762,187.4174897103089,-4.317884714201286,-161.81300927620237,-15.67934
5960738772,51.64259235864097,34.27414014713948,85.51702774336465,-37.1243587
1470173,-177.98415019794658,16.98339066468509,177.43405229412292,19.50393952
1116326},{0.,322.872632447701,4.,-154.48510511279565,-7.655531727995272,-175
.03942948598953,-1.0111563678108242,321.98759337281865,16.60122006528058,-13
3.08385014981712,-18.17273422029246,-194.634138709607,-4.033538586739716,319
.3373281801208,29.944998291495985,-110.95299226741214,-27.37121776748681,-21
3.1618091027223,-9.033977175382377,314.93636636525736},{-322.872632447701,0.
,154.48510511279565,4.,175.03942948598953,-7.655531727995272,-321.9875933728
1865,-1.0111563678108242,133.08385014981712,16.60122006528058,194.6341387096
07,-18.17273422029246,-319.3373281801208,-4.033538586739716,110.952992267412
14,29.944998291495985,213.1618091027223,-27.37121776748681,-314.936366365257
36,-9.033977175382377},{0.,688.6670919150857,4.,-114.72119024973249,-2.66535
03057498296,-650.4455719765721,-10.667982702682368,331.4293532674049,10.0696
8577122105,540.0236559683539,13.586339913460915,-511.348380335829,-20.536381
943186196,-369.6583357342561,-11.038447478035566,634.5070001073631,31.584889
64242715,158.26040739218521,2.3538027804037114,-687.2344282612197},{-688.667
0919150857,0.,114.72119024973249,4.,650.4455719765721,-2.6653503057498296,-3
31.4293532674049,-10.667982702682368,-540.0236559683539,10.06968577122105,51
1.348380335829,13.586339913460915,369.6583357342561,-20.536381943186196,-634
.5070001073631,-11.038447478035566,-158.26040739218521,31.58488964242715,687
.2344282612197,2.3538027804037114},{0.,1955.742698127284,4.,319.116359538738
95,2.610702193856547,-1851.6029756343717,-10.722043885311114,-923.3643533687
346,-9.886748303381077,1550.2743077507857,13.837048847537616,1429.2774204092
866,20.248784522486616,-1083.847095614267,-11.662601235917451,-1782.97766919
99382,-31.348033453698886,501.9941315420778,3.4859696207965745,1946.79731676
494},{-1955.742698127284,0.,-319.11635953873895,4.,1851.6029756343717,2.6107
02193856547,923.3643533687346,-10.722043885311114,-1550.2743077507857,-9.886
748303381077,-1429.2774204092866,13.837048847537616,1083.847095614267,20.248
784522486616,1782.9776691999382,-11.662601235917451,-501.9941315420778,-31.3
48033453698886,-1946.79731676494,3.4859696207965745},{0.,8267.868082849189,4
.,3927.69394068881,7.6008836161019895,-4536.125017493609,-1.1674815477135652
,-8237.513731446,-16.680749126890216,-3290.4228773517743,-17.86489003655355,
5111.248564132217,4.652781192989794,8146.673560968948,30.168278100627578,262
8.9911352802987,26.554230419975788,-5648.841605090274,-10.404652778807701,-7
996.014586041039},{-8267.868082849189,0.,-3927.69394068881,4.,4536.125017493
609,7.6008836161019895,8237.513731446,-1.1674815477135652,3290.4228773517743
,-16.680749126890216,-5111.248564132217,-17.86489003655355,-8146.67356096894
8,4.652781192989794,-2628.9911352802987,30.168278100627578,5648.841605090274
,26.554230419975788,7996.014586041039,-10.404652778807701},{0.,66714.5128785
9093,4.,49053.640102814905,11.764430373243481,5421.5038143728425,13.95034162
6292579,-41081.027085007,3.8241097110929165,-65833.36416473465,-16.221120559
06104,-55730.47653404341,-34.361001075598665,-16121.29969181952,-36.24174857
068484,32023.238065218116,-15.09440727983296,63213.194009826875,21.624753701
270276,60935.1643846494},{-66714.51287859093,0.,-49053.640102814905,4.,-5421
.5038143728425,11.764430373243481,41081.027085007,13.950341626292579,65833.3
6416473465,3.8241097110929165,55730.47653404341,-16.22112055906104,16121.299
69181952,-34.361001075598665,-32023.238065218116,-36.24174857068484,-63213.1
94009826875,-15.09440727983296,-60935.1643846494,21.624753701270276},{0.,2.1
09193656453301*^6,4.,1.931251213300796*^6,14.650157570060415,1.4274481663250
76*^6,28.242584405174718,682791.3566687256,39.659403480135936,-177073.045991
05718,43.712855697405764,-1.0070598598161684*^6,36.57092835044224,-1.6671251
576076704*^6,16.935107017831488,-2.0458959211793535*^6,-13.318096903865762,-
2.0794620445450086*^6,-49.21136470378096,-1.7621599055136922*^6},{-2.1091936
56453301*^6,0.,-1.931251213300796*^6,4.,-1.427448166325076*^6,14.65015757006
0415,-682791.3566687256,28.242584405174718,177073.04599105718,39.65940348013
5936,1.0070598598161684*^6,43.712855697405764,1.6671251576076704*^6,36.57092
835044224,2.0458959211793535*^6,16.935107017831488,2.0794620445450086*^6,-13
.318096903865762,1.7621599055136922*^6,-49.21136470378096},{0.,3.21113804748
0835*^10,4.,3.2001703830231155*^10,15.945351888106718,3.1673423099192696*^10
,35.67267128166534,3.1128780771958237*^10,62.91151228627201,3.03714973078828
5*^10,97.28714180413428,2.94067457209015*^10,138.32438810885077,2.8241116242
59851*^10,185.4525285933514,2.6882571304243427*^10,238.01117208726953,2.5340
39114531167*^10,295.2570820108493,2.362511042003884*^10},{-3.211138047480835
*^10,0.,-3.2001703830231155*^10,4.,-3.1673423099192696*^10,15.94535188810671
8,-3.1128780771958237*^10,35.67267128166534,-3.037149730788285*^10,62.911512
28627201,-2.94067457209015*^10,97.28714180413428,-2.824111624259851*^10,138.
32438810885077,-2.6882571304243427*^10,185.4525285933514,-2.534039114531167*
^10,238.01117208726953,-2.362511042003884*^10,295.2570820108493}};

b={1.4715177646857693,0.,1.7272858774111686,0.,2.4384494832359858,0.,3.86819
71196755,0.,6.1130185623361015,0.,7.4508820294899945,0.,3.393321220405215,0.
,0.03385691675272283,0.,3.3743307736651725*^-14,0.,1.23885112322841128107475
3`10.061735978331926*^-496,0.};



  • Prev by Date: Re: Returning a local array
  • Next by Date: Interesting problem: Use of NonlinearRegress inside a package
  • Previous by thread: Re: Conditionals -- what is "fastest" way to evaluate
  • Next by thread: Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]