Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 to Mathematica to isolate (or detect) that result (or badly conditioned matrix)? For instance, using the If statement, which condition should I put into the argument of If statement when my matrix is badly conditioned?

P.S. Some advice to tackle with these kind of matrices would be good.

share|improve this question
3  
You can use LinearAlgebra`MatrixConditionNumber 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 '13 at 19:28
4  
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. –  Mark McClure Feb 16 '13 at 19:36
1  
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. –  whuber Feb 16 '13 at 19:44

2 Answers 2

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, usually the choice of condition number 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 @@ Through[{Max, Min} SingularValueList[mat, Tolerance -> 0]]

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 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 will then 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 as well or 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 a matrix 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 best to stop here.

share|improve this answer

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

Quiet[Check[Inverse[N[HilbertMatrix[20]]], "Badly conditioned", Inverse::luc],   
      Inverse::luc]
share|improve this answer
    
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] –  simsurace Sep 25 '13 at 13:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.