# Solving a linear system with a badly conditioned matrix

My question is related to this one: Badly conditioned matrix (General::luc)

I'm trying to solve a linear system of the form A.x=b. As suggested in the referenced post since I'm not really interested in the inverse MMA might preform better if solve the system using x=LinearSolve[A,b], so that's what I do.

I still get a warning regarding the condition number, and the results I get later on in the computation (which are based on the solution to the linear system), do indeed show problems.

In the referenced post J.M.'s answer provides a way to determine if my matrix is indeed ill conditioned, but I don't know how can I get a better solution in case it does.

Some technical details:

1. I'm solving this equation over and over again, each time with slightly different matrix A, in order to get down the line the value for each data point, thus I need some general solution as I cannot treat each ill conditioned case separately.
2. I'm a novice on the subject, but I only guess that if I could force MMA to work with better internal accuracy I would eliminate the problem, but there's no such option as tolerance for the LinearSolve function, as there's for other linear algebra functions. Am I right? Is there a workaround?
3. A technical detail that might help: my vector b is the same for all data points, and has the form: b={0,0,0,...,0,1}, i.e. all of its entries are zero except for the last one which is unity. This means that if I were to solve the system using inversion, I wouldn't really need to find Inverse[A], but only its rightmost column. Would finding only this column be easier. I guess Cramer's rule isn't very helpful as calculating the determinant is just as bad (right?), but maybe there's a different way which I'm unaware of?
4. I don't know if it matters, but my matrix isn't that large, and has dimensions of 16X16

I shall give more details, as perhaps I wan't clear enough in the beginning:

1. I iteratively build matrices of dimensions 16X16, which hold numerical values. They all share the following structure:

• They are rather sparse.
• Most rows contain small values (with respect to the bottom one).
• The entries on the diagonal are all non-zero negative values.
• All other non-zero entries are positive.
• The bottom row is all ones.
• An example for one of my badly conditioned matrices (sorry about the formatting):
• A={{-0.04000323497710545, 0.019998382511436725, 0.019998382511436725, 1.0548487421137532*^-14, 1.0548487421137532*^-14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, {0.01999537041672573, -0.019998382511436725, 0., 0., 0., 0.03999375292816244, 2.714439127491079*^-6, 3.123535918772422*^-6, 6.247071837545091*^-6, 4.09096791281486*^-7, 9.629649721936179*^-33, 0., 0., 0., 0., 0.}, {0.01999537041672573, 0., -0.019998382511436725, 0., 0., 0.03999375292816244, 6.51319418358286*^-6, 3.123535918772422*^-6, 0., 2.857413572734654*^-6, 9.629649721936179*^-33, 0., 0., 0., 0., 0.}, {6.247071826996701*^-6, 0., 0., -1.0548487421137532*^-14, 0., 3.4666738998970245*^-33, 0.017377838870198482, 0.01999687646408122, 0.03999375292816246, 0.002619037593882741, 6.247071837545974*^-6, 0., 0., 0., 0., 0.}, {6.247071826996701*^-6, 0., 0., 0., -1.0548487421137532*^-14, 3.4666738998970245*^-33, 0.041697468145927896, 0.01999687646408122, 0., 0.01829316124631577, 6.247071837545974*^-6, 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., -0.07998750585632489, 0., 0., 0., 0., 0., 6.247071837544959*^-6, 6.247071837544959*^-6, 1.9683065157343793*^-32, 1.0839021819949165*^-32, 0.}, {0., 0., 0., 0., 0., 0., -0.059084534649437456, 0., 0., 0., 0., 0.0035343599700000386, 0.017377838870198486, 2.7144391274860995*^-6, 5.520712365254425*^-7, 0.}, {0., 0., 0., 0., 0., 0., 0., -0.03999999999999999, 0., 0., 0., 0.019996876464081232, 0.01999687646408122, 3.1235359187790123*^-6, 3.1235359187717977*^-6, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., -0.04, 0., 0., 0., 0.03999375292816246, 6.2470718375337105*^-6, 0., 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., -0.020915465350562525, 0., 0.05645626942224363, 0.0026190375938827423, 4.0909679128073285*^-7, 8.818536519796756*^-6, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.000012494143675091948, 1.698670210949542*^-31, 5.827478825726898*^-30, 0.03999375292816246, 0.03999375292816246, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.07999375292816245, 0., 0., 0., 6.247071837548034*^-6}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.07999375292816244, 0., 0., 6.247071837533809*^-6}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.04000624707183754, 0., 0.03999375292816247}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.04000624707183755, 0.03999375292816246}, {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}}
2. For each of my created matrices I want to solve the equation: $Ax=b$, where $b$ is the constant vector (always the same) b={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1} (All zeros with one at the end).

3. Currently I use LinearSolve[A,b]. As said before, for each iteration A is slightly different. For some iterations this works well, and for others (like the example I gave here), I get a warning regarding the condition number, and possible numerical error.

4. The question is as before, what can I do to improve the solutions I get. I emphasize that I the solution must be general enough that it can be incorporated into my iterations, when a warning is returned (This can by done using Check for instance to catch warnings).

# Update 2

1. Perhaps I should've written it to begin with, but the singular nature of my matrix comes as no surprise, as since one of the eigenvalues is very close to zero (perhaps numerically it is), is the reason for the ill-conditioned nature of the matrix.
2. Thanks to Daniel's and Bill's suggestions, I'm still able to get satisfactory results.
3. I'm aware of the fact that my question wasn't formulated in the best way, but I think the the solutions provided are worth keeping for others who might encounter similar issues. Thus I suggest removing the OnHold tag. How would you like me to rewrite the question?
– zhk
Apr 26, 2017 at 8:58
• What? What do mean by where? Apr 26, 2017 at 10:08
• I'm trying to solve a linear system  Where is this system you are solving?
– zhk
Apr 26, 2017 at 10:28
• There is no actual code here for anyone to attempt to improve on whatever might be the problem. Voting to close as not having sufficient data to address the question. Apr 26, 2017 at 15:52
• Something to try...since you are only making small modifications to the original A matrix, and the one you have is not burping up a warning on condition number, invert it once to get $A^{-1}$ and then use that to precondition the other runs, so if B is the next matrix and is "close" to A, solve $A^{-1} B x = A^{-1}b$ Apr 27, 2017 at 13:33

In the example (finally) provided, the matrix is rank deficient. In general not much can be done unless one knows the vector on the right is in the range space of the matrix. We'll proceed with that in mind.

mat = {{-0.04000323497710545, 0.019998382511436725,
0.019998382511436725, 1.0548487421137532*^-14,
1.0548487421137532*^-14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0.}, {0.01999537041672573, -0.019998382511436725, 0., 0., 0.,
0.03999375292816244, 2.714439127491079*^-6, 3.123535918772422*^-6,
6.247071837545091*^-6, 4.09096791281486*^-7,
9.629649721936179*^-33, 0., 0., 0., 0., 0.}, {0.01999537041672573,
0., -0.019998382511436725, 0., 0., 0.03999375292816244,
6.51319418358286*^-6, 3.123535918772422*^-6, 0.,
2.857413572734654*^-6, 9.629649721936179*^-33, 0., 0., 0., 0.,
0.}, {6.247071826996701*^-6, 0., 0., -1.0548487421137532*^-14, 0.,
3.4666738998970245*^-33, 0.017377838870198482,
0.01999687646408122, 0.03999375292816246, 0.002619037593882741,
6.247071837545974*^-6, 0., 0., 0., 0.,
0.}, {6.247071826996701*^-6, 0., 0., 0., -1.0548487421137532*^-14,
3.4666738998970245*^-33, 0.041697468145927896,
0.01999687646408122, 0., 0.01829316124631577,
6.247071837545974*^-6, 0., 0., 0., 0., 0.}, {0., 0., 0., 0.,
0., -0.07998750585632489, 0., 0., 0., 0., 0.,
6.247071837544959*^-6, 6.247071837544959*^-6,
1.9683065157343793*^-32, 1.0839021819949165*^-32, 0.}, {0., 0.,
0., 0., 0., 0., -0.059084534649437456, 0., 0., 0., 0.,
0.0035343599700000386, 0.017377838870198486,
2.7144391274860995*^-6, 5.520712365254425*^-7, 0.}, {0., 0., 0.,
0., 0., 0., 0., -0.03999999999999999, 0., 0., 0.,
0.019996876464081232, 0.01999687646408122, 3.1235359187790123*^-6,
3.1235359187717977*^-6, 0.}, {0., 0., 0., 0., 0., 0., 0.,
0., -0.04, 0., 0., 0., 0.03999375292816246,
6.2470718375337105*^-6, 0., 0.}, {0., 0., 0., 0., 0., 0., 0., 0.,
0., -0.020915465350562525, 0., 0.05645626942224363,
0.0026190375938827423, 4.0909679128073285*^-7,
8.818536519796756*^-6, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., -0.000012494143675091948, 1.698670210949542*^-31,
5.827478825726898*^-30, 0.03999375292816246, 0.03999375292816246,
0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., -0.07999375292816245, 0., 0., 0., 6.247071837548034*^-6}, {0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.07999375292816244,
0., 0., 6.247071837533809*^-6}, {0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., -0.04000624707183754, 0.,
0.03999375292816247}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., -0.04000624707183755, 0.03999375292816246}, {1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};
vec = UnitVector[Length[mat], Length[mat]];


A plausible first step in ill-conditioned cases is to work with the singular values decomposition.

{uu, ww, vv} = SingularValueDecomposition[mat];


We'll see how bad this matrix might be.

Diagonal[ww]

(* Out[325]= {4.00034774077, 0.103393916064, 0.0985619677111, \
0.0951305544638, 0.0796520702375, 0.0669929060285, 0.0548596279493, \
0.0503885069044, 0.0413276293297, 0.0400062470692, 0.0325805230661, \
0.0205858677042, 0.0199983819062, 3.61435763556*10^-6,
2.55039684488*10^-6, 0.} *)


One singular value is 0 (so it is rank deficient) and two others are small. But we might be able to simplify the linear algebra problem using this decomposition. Recalling that, up to numeric fuzz, we have uu.ww.Transpose[vv] == mat and moreover the left and right singular vector matrices are unitary (so inverse=transpose), we will transform our problem space to find x for whichww.Transpose[vv].x == Transpose[uu].vec]. First let y=Transpose[vv].x and solve for y. then we set x=vv.y.

y = LinearSolve[ww, Transpose[uu].vec];
x = vv.y

(* Out[329]= {8.44718017756*10^-10, 8.44622899399*10^-10,
8.44625619445*10^-10, 0.499999998767, 0.499999998702, \
-8.32682709879*10^-14, 8.59341503419*10^-15, 4.33994857987*10^-15,
4.6636787291*10^-15, 3.05599464256*10^-15, -3.16119352917*10^-12,
1.77271423994*10^-14, 1.04729755568*10^-14, 5.42562549792*10^-15,
5.36238940621*10^-15, 1.98712574173*10^-15} *)


We'll check this.

mat.x - vec

(* Out[332]= {1.3331704095*10^-15, -3.97230076355*10^-15, \
-4.02668630013*10^-15, 4.13664406954*10^-16, 4.84035098563*10^-16,
6.66059748146*10^-15, -2.63068456198*10^-16, 3.90350028599*10^-16,
2.32340341948*10^-16, 9.64369402065*10^-16,
4.709496088*10^-16, -1.4180482355*10^-15, -8.37760205392*10^-16, \
-1.37586298237*10^-16, -1.35056459528*10^-16, 4.4408920985*10^-16} *)


The residual is numeric fizz so we have a sound result.

One can use Chop on the result if one believes (seemingly appropriately, in this case) that the solution should be sparse.

xC = Chop[x, 10^(-8)]

(* Out[333]= {0, 0, 0, 0.499999998767, 0.499999998702, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)

• Thanks for the detailed answer. when I write the matrix symbolically it isn't singular, but I expect it to be close to singular when putting numeric values. In fact the numeric values I use are temperature dependent, and I expect the matrix to become singular for low enough temperatures. The SVD does indeed work well for me. Thanks again! May 3, 2017 at 10:14
• When does this help? For my case, I am getting "LinearSolve::nosol: Linear equation encountered that has no solution." when I solve for y. Which is the same thing I get for the original problem. ww has the same (approximate) matrix rank as mat so how does this solve the "numerical" redundancy in the matrix. Jun 8, 2022 at 13:41
• It works when the rhs is in the range space of the matrix. Jun 8, 2022 at 15:33

Turn my comment into an answer...

If you have an $A$ matrix where it works well, and the other matrix $B$ is "close" to $A$ but ill conditioned, use $A$ to precondition the problem. So solve...

$$(A^{-1} B) x = A^{-1} b$$

You only need to compute $A^{-1}$ and $A^{-1} b$ once.

As Daniel says, your A matrix is rank deficient:

Dimensions[A]
{16, 16}

Rank[A]
15


One way to proceed is to use the PseudoInverse:

s = PseudoInverse[A].b


This returns an answer (with no warnings or errors) though it will not be exact. You can test how good/bad the answer is directly:

Norm[A.s - b]
9.10073*10^-15


If this value is small enough then you can probably trust your subsequent calculations. If this grows, then you know to take some other action.

• Thanks, I added the PseudoInverse in case Check throws an warning and it works like a charm, I also check the solution as you suggested before returning the value, and so far so good. May 3, 2017 at 10:10
• I'd accept this answer too, but sadly I can accept only one of them, and Daniel's was more thorough. May 3, 2017 at 10:15

there.I have faced similar problem, My friend suggest me an approach. Using function Rationalize, Inverse and Chop.

mat = {{-0.04000323497710545, 0.019998382511436725, 0.019998382511436725, 1.0548487421137532*^-14, 1.0548487421137532*^-14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{0.01999537041672573, -0.019998382511436725, 0., 0., 0., 0.03999375292816244, 2.714439127491079*^-6, 3.123535918772422*^-6, 6.247071837545091*^-6, 4.09096791281486*^-7, 9.629649721936179*^-33, 0.,
0., 0., 0., 0.}, {0.01999537041672573, 0., -0.019998382511436725, 0., 0., 0.03999375292816244, 6.51319418358286*^-6, 3.123535918772422*^-6, 0., 2.857413572734654*^-6, 9.629649721936179*^-33, 0.,
0., 0., 0., 0.}, {6.247071826996701*^-6, 0., 0., -1.0548487421137532*^-14, 0., 3.4666738998970245*^-33, 0.017377838870198482, 0.01999687646408122, 0.03999375292816246, 0.002619037593882741,
6.247071837545974*^-6, 0., 0., 0., 0., 0.}, {6.247071826996701*^-6, 0., 0., 0., -1.0548487421137532*^-14, 3.4666738998970245*^-33, 0.041697468145927896, 0.01999687646408122, 0.,
0.01829316124631577, 6.247071837545974*^-6, 0., 0., 0., 0., 0.}, {0., 0., 0., 0., 0., -0.07998750585632489, 0., 0., 0., 0., 0., 6.247071837544959*^-6, 6.247071837544959*^-6,
1.9683065157343793*^-32, 1.0839021819949165*^-32, 0.}, {0., 0., 0., 0., 0., 0., -0.059084534649437456, 0., 0., 0., 0., 0.0035343599700000386, 0.017377838870198486, 2.7144391274860995*^-6,
5.520712365254425*^-7, 0.}, {0., 0., 0., 0., 0., 0., 0., -0.03999999999999999, 0., 0., 0., 0.019996876464081232, 0.01999687646408122, 3.1235359187790123*^-6, 3.1235359187717977*^-6, 0.},
{0., 0., 0., 0., 0., 0., 0., 0., -0.04, 0., 0., 0., 0.03999375292816246, 6.2470718375337105*^-6, 0., 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., -0.020915465350562525, 0., 0.05645626942224363,
0.0026190375938827423, 4.0909679128073285*^-7, 8.818536519796756*^-6, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.000012494143675091948, 1.698670210949542*^-31, 5.827478825726898*^-30,
0.03999375292816246, 0.03999375292816246, 0.}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.07999375292816245, 0., 0., 0., 6.247071837548034*^-6},
{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.07999375292816244, 0., 0., 6.247071837533809*^-6}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.04000624707183754, 0.,
0.03999375292816247}, {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -0.04000624707183755, 0.03999375292816246}, {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}};

vec = UnitVector[Length[mat], Length[mat]];

mat = Rationalize[mat, 10^(-30)];
xvec = N[Inverse[mat] . vec]
Chop[xvec, 10^(-9)]


ouput is

{0, 0, 0, 0.4999999987337153, 0.4999999987337153, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}


This approach worked for me. Hoping this will be helpful.

• The main issue with this method is that it hides but does not avoid ill-conditioning. A second issue, for large matrices, is that it will be slow. Sep 15, 2020 at 13:26
• @ Daniel Lichtblau, I agree with you. I thought since it worked for my case, and could be of help to others, I posted here. So far for matrix of size 20x20 for full rank matrix (which i tried) worked well. If you don't mind, Would you explain what you mean by "hides but does not avoid ill-conditioning". Sep 15, 2020 at 13:51
• If you slightly change a numerically ill-conditioned matrix by rationalizing then the exact counterpart will, in general, have full rank. So it will be invertable. You can (and do) use the inverse to solve the system. But now there is no indication of how bad that solution might be in terms of the original matrix. Using methods based on the singular values decomposition allow for assessing the ill conditioning, estimating possible error, etc. Sep 15, 2020 at 14:09
• @ Daniel Lichtblau, Thanks for the explanation. Sep 15, 2020 at 14:14