LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
Andrew Moylan
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.};
Andrew Moylan
definitions for the matrix m and vector b into Mathematica. Sorry
about that. I've put the definitions of m and b into a Mathematica
notebook that you can get from http://andrew.j.moylan.googlepages.com/files.
This might be more convenient.
Andrew Moylan
On Dec 18, 11:18 am, "Andrew Moylan" <andrew.j.moy...@gmail.com>
wrote:
mcmc...@unca.edu
wrote:
> 1. LinearSolve[m, b] is not equivalent to LinearSolve[m][b].
Well, I was definitely hoping someone might have clarified
this interesting issue. I haven't totally figured it out
myself, but I've got a few observations. My basic guess is
that LinearSolve[m][b] introduces the possibility of
compounded error via intermediate matrix factorizations that
a direct LinearSolve[m,b] avoids. Here's a simple example
illustrating this.
Needs["ComputerArithmetic`"];
m3 = Map[ComputerNumber, HilbertMatrix[3], {2}];
b3 = Table[ComputerNumber[1], {3}];
ans1 = LinearSolve[m3, b3, Method -> "Direct"];
ans2 = LinearSolve[m3, Method -> "Direct"][b3];
Max[Abs[ans1 - ans2]]
0.1700000000000000000
This example might seem contrived, in that we've used a low
precision approximation to an ill-conditioned matrix. It
does illustrate the basic point that LinearSolve[m][b] and
LinearSolve[m,b] are different. Also, it's interesting to
note that ans2 is relatively close to the following:
{lu, p, c} = LUDecomposition[m3];
l = lu SparseArray[{i_, j_} /; j < i -> 1, {3, 3}] +
IdentityMatrix[3];
u = lu SparseArray[{i_, j_} /; j >= i -> 1, {3, 3}];
y = LinearSolve[l, b3[[p]],
Method -> "Direct"];
ans3 = LinearSolve[u, y,
Method -> "Direct"];
Max[Abs[ans3 - ans2]]
0.01000000000000000000
If I we're going to implement something like a Direct
LinearSolveFunction from scratch, I suppose I'd use an
LUDecompostion. Perhaps the error can arise here?
In fact, the LUDecomposition is part of the
LinearSolveFunction returned by LinearSolve, when the
Method is "Direct":
lu1 = LinearSolve[m3, Method -> "Direct"][[2, 3]];
lu2 = LUDecomposition[m3];
lu1 === lu2
Now, in your example, it appears that the method is
"IterativeRefinement". I figured this out simply by playing
with the possible methods and looking at the answers.
Here's an easy way to generate a list of possible Methods:
LinearSolve[m, b, Method -> bogus];
Strangely, the "IterativeRefinement" method is not
documented in the LinearSolve documentation or the
MatrixComputations documentation; rather, it is a
LeastSquares method.
Of course, this is all conjecture; I have no direct
knowledge of the LinearSolve algorithms.
Mark
Daniel Lichtblau
>
>
> this interesting issue. I haven't totally figured it out
> myself, but I've got a few observations. My basic guess is
> that LinearSolve[m][b] introduces the possibility of
> compounded error via intermediate matrix factorizations that
> a direct LinearSolve[m,b] avoids. Here's a simple example
> illustrating this.
Strangely enough, the issue is that LinearSolve[m][b] gives a good
result, whereas LinearSolve[m, b] gives garbage. I'm fairly certain this
is simply a bug. It was filed as such, and will be investigated.
Daniel Lichtblau
Wolfram Research