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

• To: mathgroup at smc.vnet.net
• Subject: [mg84356] Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
• From: mcmcclur at unca.edu
• Date: Thu, 20 Dec 2007 00:05:46 -0500 (EST)
• References: <fk73ks\$65f\$1@smc.vnet.net>

On Dec 17, 7:18 pm, "Andrew Moylan" <andrew.j.moy... at gmail.com>
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, {2}];
b3 = Table[ComputerNumber, {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;
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

• Prev by Date: Re: FileNameSetter button return state
• Next by Date: Re: FileNameSetter button return state
• Previous by thread: Re: LinearSolve[m, b] is not equivalent to LinearSolve[m][b]
• Next by thread: Re: Re: LinearSolve[m, b] is not equivalent to LinearSolve