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].