Eigenvalues and Determinant of a large matrix

Question

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?

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.