Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

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

Search the Archive

LinearSolve, precision and matrix conditioning

  • To: mathgroup at smc.vnet.net
  • Subject: [mg87312] LinearSolve, precision and matrix conditioning
  • From: "Andrew Moylan" <andrew.j.moylan at gmail.com>
  • Date: Mon, 7 Apr 2008 05:18:10 -0400 (EDT)

>From tutorial/LinearAlgebraMatrixComputations under "Estimating and
Calculating Accuracy", an example of applying LinearSolve to an
ill-conditioned matrix at machine precision:

mat = N[{{1, 2, 3}, {4, 5, 6},
    {3, 6, 9 + 10^11 $MachineEpsilon}}];
mat2 = mat.mat;
vec = {3, 4, 5};
sol = LinearSolve[mat2, vec];
Norm[mat2.sol - vec]
>> 0.000304962

According to tutorial/LinearAlgebraMatrixComputations, a better result is
obtained using arbitrary precision arithmetic:

matAcc = SetPrecision[mat2, 16];
sol = LinearSolve[matAcc, vec];
Norm[matAcc.sol - vec]
>> 0

1. I don't think this shows anything, because that 0 arose from a
low-precision subtraction:

matAcc.sol // Precision
>> 3.63143

Do you agree?

2. Nonetheless, by increasing the precision of matAcc.sol, I think we find
that the result obtained using arbitrary precision input (at precision 16)
*is* actually better:

Norm[SetPrecision[matAcc, 50].SetPrecision[sol, 50] - vec]
>> 1.759196429592308379*10^-18

(Increasing the precision of mat2.sol in the machine precision example above
does *not* improve the Norm of the error.)

I think it is strange that the arbitrary precision LinearSolve (at precision
16) should produce a much more accurate result than the machine precision
LinearSolve. In fact, if we set matAcc *lower* than machine precision (e.g.,
matAcc = SetPrecision[mat2, 14]) LinearSolve still appears to produce a
superior result.

Any ideas on why?

3. A previous curious dependence of LinearSolve on its input precision:

http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/
thread/8a58b79fb6c80d23/7d6af6f6e084d0e8

I don't think this issue is quite the same, because there does not seem to
be any difference when using LinearSolve[m][b] instead of LinearSolve[m, b].



  • Prev by Date: Re: CreateDialog Pecularity
  • Next by Date: Re: A small syntax error causing Kernel crash, Version 6.0.1
  • Previous by thread: Re: Printing Mathematica stuff via KPrinter
  • Next by thread: Re: Poisson equation with boundary conditions on rectangle