Services & ResourcesWolfram ForumsMathGroup Archive

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

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[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

• Follow-Ups:
  • Re: Re: LinearSolve[m, b] is not equivalent to LinearSolve
    • From: Daniel Lichtblau <danl@wolfram.com>