MathGroup Archive 2007

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

Search the Archive

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

  • To: mathgroup at smc.vnet.net
  • Subject: [mg84315] Re: [mg84283] LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
  • From: DrMajorBob <drmajorbob at bigfoot.com>
  • Date: Wed, 19 Dec 2007 04:08:47 -0500 (EST)
  • References: <17656120.1197990404317.JavaMail.root@m35>
  • Reply-to: drmajorbob at bigfoot.com

I get different answers than you, probably because your m and b were not 
PRECISELY recorded in your e-mail.

sol1 = LinearSolve[m, b]
m.sol1 - b // Norm

{-4.38029*10^45, -4.67816*10^43, 1.69885*10^54, -4.52491*10^54,
  5.72302*10^49, -1.85773*10^55, -6.2527*10^52, 2.94226*10^53,
  1.3215*10^53, -1.43897*10^52, -6.97568*10^54, 2.42227*10^52,
  5.90221*10^53, -3.03106*10^48, 2.63164*10^54, -9.48003*10^51,
  9.50428*10^51, 7.73519*10^44, 2.0085*10^54, -7.66248*10^53}

1.76753*10^64

sol2 = LinearSolve[m][b]
m.sol2 - b // Norm

LinearSolve::luc: Result for LinearSolve of badly conditioned matrix \
{{0.,120.,4.,-120.,-16.,120.,36.,-120.,-64.,120.,<<10>>},<<9>>,<<10>>}\
  may contain significant numerical errors. >>

{-0.000039224, -1.15304*10^-7, 47613.2, -1.47015*10^6, 2.05617, \
-1.46211*10^6, -1687.25, 158329., 4747.89, -1132.38, -238101., \
13081.9, 21205.6, -0.169042, 85356.7, -719.205, 341.471, \
0.0000905787, 73034.3, -53218.9}

400.012

Here are another wild/tame pair:

n = Rationalize[m, 0];
LinearSolve[n, b]
n.% - b // Norm

{-1.19108*10^45, -1.8537*10^44, 9.46945*10^53, -4.73058*10^55,
  4.98366*10^49, -2.71606*10^55, -3.63477*10^52, -1.25435*10^55,
  1.15077*10^53, -2.10438*10^52, -5.81652*10^54, -1.03688*10^54,
  5.13971*10^53, -6.87898*10^48, 2.03201*10^54, 1.46442*10^52,
  8.27641*10^51, 2.81116*10^45, 1.61668*10^54, 1.5325*10^54}

1.30003*10^64

LinearSolve[n][b]
n.% - b // Norm

{-0.0331528, 0.340685, 47617.5, -1.47014*10^6, 2.05623, -1.4621*10^6, \
-1687.3, 158330., 4748.02, -1132.38, -238111., 13081.9, 21206.2, \
-0.169041, 85360., -719.203, 341.481, 0.0000905786, 73036.2, \
-53218.5}

6.5775*10^18

All four errors are miniscule when compared to the determinant:

Det@m

-1.14512*10^187

Bobby

On Mon, 17 Dec 2007 18:16:03 -0600, Andrew Moylan  
<andrew.j.moylan at gmail.com> wrote:

> 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,13=
6.9873
> 54117796,8.062099516267708,-120.51022451941246,-44.15773166836576},{0.=
,190.5
> 9730297009568,4.,-140.79278017868444,-11.819078485136764,17.4078119139=
64213,
> 14.191990544591762,115.0747432704792,-4.317884714201286,-187.417489710=
3089,-
> 15.679345960738772,161.81300927620237,34.27414014713948,-51.6425923586=
4097,-
> 37.12435871470173,-85.51702774336465,16.98339066468509,177.98415019794=
658,19
> .503939521116326,-177.43405229412292},{-190.59730297009568,0.,140.7927=
801786
> 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.1=
243587
> 1470173,-177.98415019794658,16.98339066468509,177.43405229412292,19.50=
393952
> 1116326},{0.,322.872632447701,4.,-154.48510511279565,-7.65553172799527=
2,-175
> .03942948598953,-1.0111563678108242,321.98759337281865,16.601220065280=
58,-13
> 3.08385014981712,-18.17273422029246,-194.634138709607,-4.0335385867397=
16,319
> .3373281801208,29.944998291495985,-110.95299226741214,-27.371217767486=
81,-21
> 3.1618091027223,-9.033977175382377,314.93636636525736},{-322.872632447=
701,0.
> ,154.48510511279565,4.,175.03942948598953,-7.655531727995272,-321.9875=
933728
> 1865,-1.0111563678108242,133.08385014981712,16.60122006528058,194.6341=
387096
> 07,-18.17273422029246,-319.3373281801208,-4.033538586739716,110.952992=
267412
> 14,29.944998291495985,213.1618091027223,-27.37121776748681,-314.936366=
365257
> 36,-9.033977175382377},{0.,688.6670919150857,4.,-114.72119024973249,-2=
.66535
> 03057498296,-650.4455719765721,-10.667982702682368,331.4293532674049,1=
0.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},{-6=
88.667
> 0919150857,0.,114.72119024973249,4.,650.4455719765721,-2.6653503057498=
296,-3
> 31.4293532674049,-10.667982702682368,-540.0236559683539,10.06968577122=
105,51
> 1.348380335829,13.586339913460915,369.6583357342561,-20.53638194318619=
6,-634
> .5070001073631,-11.038447478035566,-158.26040739218521,31.584889642427=
15,687
> .2344282612197,2.3538027804037114},{0.,1955.742698127284,4.,319.116359=
538738
> 95,2.610702193856547,-1851.6029756343717,-10.722043885311114,-923.3643=
533687
> 346,-9.886748303381077,1550.2743077507857,13.837048847537616,1429.2774=
204092
> 866,20.248784522486616,-1083.847095614267,-11.662601235917451,-1782.97=
766919
> 99382,-31.348033453698886,501.9941315420778,3.4859696207965745,1946.79=
731676
> 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.86808284=
9189,4
> .,3927.69394068881,7.6008836161019895,-4536.125017493609,-1.1674815477=
135652
> ,-8237.513731446,-16.680749126890216,-3290.4228773517743,-17.864890036=
55355,
> 5111.248564132217,4.652781192989794,8146.673560968948,30.1682781006275=
78,262
> 8.9911352802987,26.554230419975788,-5648.841605090274,-10.404652778807=
701,-7
> 996.014586041039},{-8267.868082849189,0.,-3927.69394068881,4.,4536.125=
017493
> 609,7.6008836161019895,8237.513731446,-1.1674815477135652,3290.4228773=
517743
> ,-16.680749126890216,-5111.248564132217,-17.86489003655355,-8146.67356=
096894
> 8,4.652781192989794,-2628.9911352802987,30.168278100627578,5648.841605=
090274
> ,26.554230419975788,7996.014586041039,-10.404652778807701},{0.,66714.5=
128785
> 9093,4.,49053.640102814905,11.764430373243481,5421.5038143728425,13.95=
034162
> 6292579,-41081.027085007,3.8241097110929165,-65833.36416473465,-16.221=
120559
> 06104,-55730.47653404341,-34.361001075598665,-16121.29969181952,-36.24=
174857
> 068484,32023.238065218116,-15.09440727983296,63213.194009826875,21.624=
753701
> 270276,60935.1643846494},{-66714.51287859093,0.,-49053.640102814905,4.=
,-5421
> .5038143728425,11.764430373243481,41081.027085007,13.950341626292579,6=
5833.3
> 6416473465,3.8241097110929165,55730.47653404341,-16.22112055906104,161=
21.299
> 69181952,-34.361001075598665,-32023.238065218116,-36.24174857068484,-6=
3213.1
> 94009826875,-15.09440727983296,-60935.1643846494,21.624753701270276},{=
0.,2.1
> 09193656453301*^6,4.,1.931251213300796*^6,14.650157570060415,1.4274481=
663250
> 76*^6,28.242584405174718,682791.3566687256,39.659403480135936,-177073.=
045991
> 05718,43.712855697405764,-1.0070598598161684*^6,36.57092835044224,-1.6=
671251
> 576076704*^6,16.935107017831488,-2.0458959211793535*^6,-13.31809690386=
5762,-
> 2.0794620445450086*^6,-49.21136470378096,-1.7621599055136922*^6},{-2.1=
091936
> 56453301*^6,0.,-1.931251213300796*^6,4.,-1.427448166325076*^6,14.65015=
757006
> 0415,-682791.3566687256,28.242584405174718,177073.04599105718,39.65940=
348013
> 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.21113=
804748
> 0835*^10,4.,3.2001703830231155*^10,15.945351888106718,3.16734230991926=
96*^10
> ,35.67267128166534,3.1128780771958237*^10,62.91151228627201,3.03714973=
078828
> 5*^10,97.28714180413428,2.94067457209015*^10,138.32438810885077,2.8241=
116242
> 59851*^10,185.4525285933514,2.6882571304243427*^10,238.01117208726953,=
2.5340
> 39114531167*^10,295.2570820108493,2.362511042003884*^10},{-3.211138047=
480835
> *^10,0.,-3.2001703830231155*^10,4.,-3.1673423099192696*^10,15.94535188=
810671
> 8,-3.1128780771958237*^10,35.67267128166534,-3.037149730788285*^10,62.=
911512
> 28627201,-2.94067457209015*^10,97.28714180413428,-2.824111624259851*^1=
0,138.
> 32438810885077,-2.6882571304243427*^10,185.4525285933514,-2.5340391145=
31167*
> ^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.393321220405=
215,0.
> ,0.03385691675272283,0.,3.3743307736651725*^-14,0.,1.23885112322841128=
107475
> 3`10.061735978331926*^-496,0.};
>
>
>



-- =

DrMajorBob at bigfoot.com


  • Prev by Date: Re: Updating NumberPadding?
  • Next by Date: Re: Mathematica SIG (Washington D.C. and Northern
  • Previous by thread: Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
  • Next by thread: Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]