       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: