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.

Can anybody kindly explain to me what is going wrong regarding a simple problem I have? I can find the eigenvalues of a large matrix using Eigenvalues[], but when I try to find the determinant of the matrix minus eigenvalue times identity matrix, the result is not equal to zero. To make it clear please take a look at the following example:

n = 150;
m = RandomReal[{-5, 5}, {n, n}];
EIG = Eigenvalues[m];
DET = Det[(m - (EIG[[n]])*IdentityMatrix[n])]

For small values of n, DET equals zero, but when n becomes larger DET is no longer zero, and for n = 150 (as in the example) DET will be a huge number!

Could you please tell me whether the problem is from Eigenvalues[] or Det[] in Mathematica, and how I can sort it out?

share|improve this question
1  
Try it with m = RandomReal[{-5, 5}, {n, n}, WorkingPrecision -> 200]; Even at that it will only give a few digits of accuracy in the result. MatrixRank will be more reliable (even at machine precision). –  Daniel Lichtblau Feb 27 '13 at 17:12
    
Conversely, symbolic calculation does work if your matrix is not too big :n = 25; m = RandomReal[{-5, 5}, {n, n}] // Rationalize[#, 10^-10] &; EIG = Eigenvalues[m]; DET = Det[(m - (EIG[[n]])*IdentityMatrix[n])] // N –  chris Feb 27 '13 at 17:29
3  
It makes sense only to evaluate this determinant relative to the determinant of $m$. When you divide DET by Det[m] you will get values around $10^{-13}$ down to $10^{-16}$ in size, indicating that only one to four decimals of accuracy have actually been lost out of the original $17$. In other words, DET is just the opposite of a "huge number": it is extremely small, relative to a natural measure of numerical size for this problem. –  whuber Feb 27 '13 at 18:54
add comment

1 Answer 1

It's just a matter of the difficulty inherent in the numerical computation of determinants. Here's what Cleve Moler has to say about determinants and characteristic polynomials in chapter 10 of his book on numerical computing:

Like the determinant itself, the characteristic polynomial is useful in theoretical considerations and hand calculations, but does not provide a sound basis for robust numerical software.

Moreover, the condition number of your matrix is, uh, huge:

n = 150;
m = RandomReal[{-5, 5}, {n, n}]; 
LinearAlgebra`MatrixConditionNumber[m]

(* Out: 9000.84 *)

Edit

You can certainly compute the determinant of an exact matrix with several hundred entries.

SeedRandom[1];
size = 300;
m = RandomInteger[25, {size, size}];
Det[m] // AbsoluteTiming

(* Out: {0.385179, 
-31619667635191239031906201038858595859003061431355557157251141946641068769797\
004846355898301331685271103794188688002223539528882633242121673402163964281634\
084344740433631930627048698258792079669124806129185893568239852670285556224752\
376916371921006623721395633460341811343837292528355079446388064508960057443123\
056664896757423781301407944339347522658235223949579536040655676063666110229006\
390862616701174010869791324107235491084936157417121971773972255891505320775553\
355811238993293311257225265815528319828336822493711170605717238245535722011317\
88870176752254707662515496}
*)

We can force the determinant to be zero, by setting two rows equal to one another to create a singular matrix. We can then compare the exact computation Det[m] to the computation with approximate numbers Det[N[m]].

m[[-1]] = m[[1]];
Det[m]
Det[N[m]]

(* Out: 
  0
  -8.642842537204658*10^553
*)

The determinant of the numerical matrix is very far off, even though the entries are floating point integers. Now, the condition number is effectively infinite, since the matrix is singular.

LinearAlgebra`MatrixConditionNumber[N[m]]

(* Out: 3.46024*10^17 *)

Even though you can compute the determinants of such matrices, my advice is still don't. Determinant computation is expensive, requiring on the order of $n^3$ computations for an $n\times n$ matrix. Furthermore, there's virtually no large-scale problem that cannot be recast without reference to the determinant.

Eigenvalue computation for exact matrices is much worse, as the determinant is just one of many coefficients in the characteristic polynomial. Numerical eigenvalue computation is typically more stable than that for the determinant simply because, even if the matrix is singular, you've probably got eigenvalues far from zero. In terms of stability, though, the eigenvalues have nothing on the singular values, which can often be used instead.

If you want to learn about these issues in detail, then there are many good books on numerical linear algebra. A fantastic, intro level text written from the numerical perspective is Matrix Analysis by Carl Meyer. A more advanced text, and bound to be a classic, is Numerical Linear Algebra by Trefethen and Bau.

share|improve this answer
1  
Caveat: It can be useful numerically at low dimension for, say, area or volume integrals, or other operations that involve figuring out a volume element. Also useful in some computational geometry, again in low dimension (2 or 3 comes to mind as common cases for comp geom...) –  Daniel Lichtblau Feb 27 '13 at 17:34
    
@DanielLichtblau Dimension = 1 works, too. –  Mark McClure Feb 27 '13 at 17:36
    
How soon I forget... –  Daniel Lichtblau Feb 27 '13 at 17:37
1  
What are you considering exact? The elements of the matrix themselves aren't exact, they are machine precision numbers. –  chuy Feb 28 '13 at 14:10
3  
To amplify slightly @chuy's point, the "exact" result must really be regarded as an interval range. The numeric instability of determinants is such that this range is huge. The eigenvalues are under better control; check standard numerical linear algebra refs to see what to expect there. –  Daniel Lichtblau Feb 28 '13 at 14:46
show 3 more comments

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.