Badly conditioned matrix (General::luc)

Question

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.

Answer 1

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][[1]], {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 LinearAlgebra`MatrixConditionNumber[]. 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[LinearAlgebra`MatrixConditionNumber[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.

Answer 2

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

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

Answer 3

Just a note on computing the condition number in Mathematica:

Besides LinearAlgebra`MatrixConditionNumber[] (in V11.2+, use LinearAlgebra`Private`MatrixConditionNumber[] 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][[2]]], 
   x_ /; x == 0];
lu[mat_] := Quiet@ Last@ LUDecomposition[mat];
BenchmarkPlot[{LinearAlgebra`MatrixConditionNumber, svl, svd, lu}, 
 N[HilbertMatrix[#]] &, {10, 30, 100, 300, 1000, 3000, 10000}]