Why the sparse array is ill-conditioned?
Let me present another workaround. By default, LinearSolve[]
is using a multifrontal method (Method -> "Multifrontal"
) on the given SparseArray[]
, and it seems it is having trouble with the internal thresholding, which results in the erroneous error message.
Thus, you might consider using a different method instead. In particular, the structure of the matrix involved suggests the use of Method -> "Banded"
:
sm = SparseArray[mat];
ls = LinearSolve[sm, Method -> "Banded"];
which gives a pretty good result in this case:
BlockRandom[SeedRandom[42];
b = sm.(x = RandomReal[1, 110]);
Norm[lsx[b] - x, ∞]/Norm[x, ∞]]
2.2219*10^-16
Here is a histogram similar to Michael's for the "Banded"
strategy, with the slight replacement of using the max-norm instead of the 2-norm:
Here is what appears to be a slight improvement on Mr. Wizard's workaround. How much to Chop[]
should probably depend on both the magnitude (norm) and precision of the matrix. Something like this:
LinearSolve[SparseArray@Chop[mat, Norm[mat] $MachineEpsilon]]
With only one test example, and given that this appears to be a bug, it's hard to test the robustness of this approach. Note that lsM = LinearSolve[mat]
is accurate to machine precision.
lsM = LinearSolve[mat];
b = mat.(x0 = RandomReal[1, 110]);
Norm[lsM@b - #]/Norm[#] &[x0]
(* 2.07362*10^-16 *)
Chopping to machine precision produces a well-conditioned matrix/linear-solve-function, whose error is only a small multiple of lsM
:
lsChop = LinearSolve[SparseArray@Chop[mat, Norm[mat] $MachineEpsilon]];
Table[(b = mat.(x0 = RandomReal[{-1, 1}, 110]);
Norm[lsChop@b - #]/Norm[#] &[x0])/(Norm[lsM@b - #]/Norm[#] &[x0]),
{5000}] //
Histogram[#, {0.2}, PlotRange -> All,
PlotLabel -> "Error ratio (SA@Chop[$MachineEpsilon]/mat)"] &
The default tolerance for Chop[]
has larger error (by a factor on the order of 1000):
lsChop2 = LinearSolve[SparseArray@Chop@mat];
Table[(b = mat.(x0 = RandomReal[{-1, 1}, 110]);
Norm[lsChop2@b - #]/Norm[#] &[x0])/(Norm[lsM@b - #]/Norm[#] &[x0]),
{5000}] //
Histogram[#, {100}, PlotRange -> All,
PlotLabel -> "Error ratio (SA@Chop[default]/mat)"] &