With some matrices, I am receiving the following message:

Inverse::luc Result for Inverse of badly conditioned matrix (M) may contain significant numerical errors.


How can I tell Mathematica to isolate (or detect) that result (or badly conditioned matrix)? For instance, using the If statement, which condition should I put into the first argument of If when my matrix is badly conditioned?

P.S. Some advice to tackle matrices like these would be good.

• You can use LinearAlgebraMatrixConditionNumber to find the condition number of your matrix. What value you choose as the threshold depends on your particular application... In general, a condition number that is $\mathcal{O}(10^n)$ can make you lose up to $n$ digits of accuracy (in addition to FP errors).
– rm -rf
Feb 16, 2013 at 19:28
• Numerical inversion of a matrix is both slow and poorly conditioned, a fact mentioned in most books on numerical analysis. Typically, a solution to an applied problem can be formulated in such a way that matrix inversion is expressed in terms of the solution of a system of equations. Naturally, it would be much easier to provide details if you provided more details of your problem. Feb 16, 2013 at 19:36
• If you already happen to know the eigenvalues or singular values of the matrix, then compare the one of smallest absolute value to the one of largest absolute value: that ratio is the reciprocal of the condition number. Wikipedia gives some other alternative calculations for special circumstances. Feb 16, 2013 at 19:44

This might be as good a time as any to distill the collective wisdom of Messrs. Huber, McClure, and Toad R. M.

As already mentioned, there is this quantity of great interest to people in the business of solving simultaneous linear equations, called the condition number, and conventionally denoted by the symbol $\kappa$. This is usually associated with a matrix $\mathbf A$ that figures as the matrix of coefficients in the linear system, and thus, we have the symbol $\kappa(\mathbf A)$.

Further, there is not just one condition number, but a number of them, all dependent on the underlying matrix norm used in their definition. One then speaks of the "2-norm condition number", $\kappa_2(\mathbf A)$, the "$\infty$-norm condition number", $\kappa_\infty(\mathbf A)$... and so on.

Now, why should we be interested in this condition number? One for instance cannot depend on the determinant, as it is not a reliable measure of how badly a coefficient matrix will behave in a linear system (see also this answer). In inexact arithmetic, the condition number is pretty much the only nice diagnostic you have.

For applications, the choice of condition number usually does not matter much, since if some matrix $\mathbf A$ has a large/small $p$-norm condition number, the $q$-norm condition number of the same matrix will be of comparable magnitude. In any event, one usually wants matrices whose condition numbers are as near to unity as can be, as this is the case of well-conditioning. Conversely, if your matrix's condition number is "huge" (for some application-dependent definition of "huge"), then your matrix might as well be singular. (A matrix that is truly singular has a condition number of $\infty$.)

The two-norm condition number is easily computed in Mathematica:

cond2[mat_?MatrixQ] :=
Divide @@ Table[SingularValueList[mat, k, Tolerance -> 0][], {k, {1, -1}}]


As you might ascertain from this routine, this condition number requires the computation of the singular value decomposition (SVD). This is quite expensive, and sometimes one instead uses condition number estimators, which require much less computational effort. (One thing it can do that the next method I am about to describe can't is that it can be applied to non-square matrices as well.)

Mathematica has a condition number estimator built-in, in the form of the undocumented function LinearAlgebraMatrixConditionNumber[]. This function can be set to estimate either the $1$-norm or $\infty$-norm condition number, depending on the setting of its Norm option. The Hager-Higham condition estimator, which is the underlying algorithm, almost always gives a result that is equal to or near the exact condition number, though there are matrices that can defeat it. (Luckily, these counterexamples are rather contrived and do not seem to crop up in practice.)

So, what to do in Mathematica? As R. M. notes, this is application-dependent, but you will want to use the heuristic that when solving a linear system with coefficient matrix $\mathbf A$, you stand to lose $\approx\log_b(\kappa(\mathbf A))$ base-$b$ digits in the solution of your linear system. You thus want to do something like If[LinearAlgebraMatrixConditionNumber[A] < 1/\$MachineEpsilon, (* code for well-conditioned case *), (* code for ill-conditioned case *)], to use a typical example.

Mark's advice is only slightly related to the matter at hand, but he is of course right: one often does not need to compute inverses, and what an inverse can do, an appropriate matrix decomposition can do just as well or even better. (There are instances where inverses are genuinely needed, like variance-covariance matrices, but they are few and far between.) For instance, if you are doing something like

x = Inverse[A].b


this is better done using the functionality of LinearSolve[], which internally stores an LU decomposition:

lf = LinearSolve[A];
x = lf[b]


There are still a number of details to say (or I have forgotten), but this post is getting too long already, and I think it's best to stop here.

• Hi, great answer! Just a comment. In the definition of cond2, you forgot the square brackets around the evaluation of SingularValueList. Nov 18, 2015 at 8:43
• Yes, will fix later. Thank you very much! Nov 18, 2015 at 8:46

Use Check with 3 arguments, possibly in conjunction with Quiet. For example:

Quiet[Check[Inverse[N[HilbertMatrix]], "Badly conditioned", Inverse::luc],
Inverse::luc]

• Inspired by Murray, you can define your own function which checks whether a matrix is badly conditioned as follows: BadlyConditionedQ[mat_?MatrixQ]:=If[Quiet[Check[Inverse[mat], True, Inverse::luc], Inverse::luc]===True,True,False] Sep 25, 2013 at 13:48

Just a note on computing the condition number in Mathematica:

Besides LinearAlgebraMatrixConditionNumber[] (in V11.2+, use LinearAlgebraPrivateMatrixConditionNumber[] instead), LUDecomposition[] returns the ∞-norm condition number as the third element of the returned list. It seems somewhat faster than MatrixConditionNumber[], plus you get the LU decomposition to use when appropriate.

Similarly, SingularValueDecomposition[] may be used instead of SingularValueList[] to compute the 2-norm condition number. It is slightly slower than SingularValueList[], but again you get the decomposition to use. It's probably not worth checking SingularValueList[] before computing SingularValueDecomposition[] in most applications where the SVD might be conditionally used.

Thus, it seems to make sense to use the normal functions such as LinearSolve[] as recommended by @J.M., and use Check and Quiet as recommended by @murray to react to ill-conditioned problems as necessary. It doesn't seem likely that one could speed up code by first checking the condition number before calling LinearSolve[], for instance. (LinearSolve[] takes about the same time as LUDecomposition[]. Indeed, judging from the LAPACK routines, the matrix has to be factored before the condition number can be computed.) Needs["GeneralUtilities"];
Clear[svl, svd, lu];
svl[mat_] := First@ Ratios@ MinMax@ SingularValueList[mat];
svd[mat_] := Ratios@ MinMax@
DeleteCases[
Diagonal[SingularValueDecomposition[mat][]],
x_ /; x == 0];
lu[mat_] := Quiet@ Last@ LUDecomposition[mat];
BenchmarkPlot[{LinearAlgebraMatrixConditionNumber, svl, svd, lu},
N[HilbertMatrix[#]] &, {10, 30, 100, 300, 1000, 3000, 10000}]

• I should add that LinearSolve[] actually computes the condition number internally; one can access it like so: ls = LinearSolve[N[HilbertMatrix]]; ls["ConditionNumber"]. Jan 29, 2017 at 17:16
• @J.M. I assume it stores the complete output of LUDecomposition, at least when appropriate. An internal Trace won't confirm that though. But thanks, I didn't know about ls["ConditionNumber"]. Jan 29, 2017 at 17:38
• Yes, LinearSolveFunction[] internally stores the packed-together lower and upper triangular matrices, which is the same format returned by LUDecomposition[]. Jan 29, 2017 at 17:40
• @J.M. And there's "getL" and "getU", too. But "getPermutations" does not seem to work. Jan 29, 2017 at 17:48