So You Need the Eigenvalues of a Matrix
Most people hit this problem when they're working with systems of differential equations or trying to diagonalize a matrix for a simulation. The direct approach is straightforward enough in theory. You set up the characteristic equation det(A - I) = 0 and solve for . That's it on paper. For a 2x2 matrix, you get a quadratic. Just use the quadratic formula and you're done. For a 3x3, it's a cubic. You can use Cardano's formula, but honestly that's more work than it's worth most of the time. If the coefficients are nice integers, you might spot a rational root by trial. If they're not, numerical methods take over. For a 4x4 and above, you stop trying to do it analytically. The degree of the characteristic polynomial hits five or higher and there's no general algebraic solution. I've seen people try to factor these by hand in grad school. Nobody comes out of that experience feeling good about it.
Here's the thing nobody tells you early on: computing the characteristic polynomial explicitly is often the wrong move. It's numerically unstable. The coefficients can blow up even for moderate-sized matrices because you're expanding determinants that involve products of many entries. A 10x10 matrix's characteristic polynomial will have coefficients that vary across dozens of orders of magnitude, and floating point arithmetic eats that alive.
The QR Algorithm Is What You Actually Use
In practice, the standard method for finding eigenvalues of a dense matrix is the QR algorithm with shifts. You iterate through a sequence of similarity transformations A_k = Q_k R_k, where Q_k is orthogonal and R_k is upper triangular, then set A_{k+1} = R_k Q_k. The diagonal entries converge to the eigenvalues. With an appropriate shift strategy, convergence is typically cubic for distinct eigenvalues. The shifts matter a lot. Without shifts, the basic QR iteration can be painfully slow, especially when eigenvalues have similar magnitudes. The implicit double shift is the industry standard because it handles complex conjugate pairs without ever leaving the real domain. LAPACK's DGEEV and ZGEEV routines implement exactly this, and they're what almost every high-level language calls under the hood. If your matrix is symmetric or Hermitian, use the divide-and-conquer or QR variant for symmetric tridiagonal matrices instead. LAPACK's DSYEV/DSYEVR routines reduce to tridiagonal form first via Householder reflections, then solve the smaller problem. This is significantly faster and more accurate than the general nonsymmetric case. The condition number of the eigenvalue problem for symmetric matrices is 1, meaning the eigenvalues are as well-conditioned as they can possibly be. Nonsymmetric matrices don't have that luxury.
Get the Full Details
A Problem I Actually Had With This
Last year I was debugging a structural mechanics code that needed the eigenvalues of a stiffness matrix. The matrix was 6x6, symmetric, and theoretically positive definite. The eigenvalues came back, but two of them were slightly negative — on the order of 1e-12 relative to the largest eigenvalue. The matrix was fine, but floating point roundoff in the Cholesky factorization inside the eigensolver introduced the corruption. The workaround was to symmetrize the input explicitly before calling the solver. I added a step that replaced the matrix with (A + A^T)/2, which eliminated the tiny antisymmetric component that was destabilizing the computation. After that, all eigenvalues were cleanly positive. It sounds trivial but I wasted about three hours chasing a phantom bug before I realized the matrix I was feeding in wasn't exactly symmetric due to accumulation error from previous operations.
When Standard Methods Break Down
There are real scenarios where the QR approach is the wrong tool. If your matrix is sparse and very large, running a dense eigensolver on it is wasteful and sometimes impossible due to memory. In that case, you'd want Lanczos or Arnoldi iteration, which only require matrix-vector products rather than full factorizations. ARPACK is the standard library here, and it computes a small number of eigenvalues — usually the largest or smallest in magnitude — rather than the full spectrum. Another case: if your matrix is defective, meaning it doesn't have a full set of linearly independent eigenvectors, the eigenvalues themselves are still computable, but any downstream use of the eigenvector basis will break. Jordan normal form exists for defective matrices, but computing it numerically is ill-conditioned and not recommended. If you encounter this, you're usually dealing with a problem that needs a different formulation altogether. Also worth noting: condition numbers for eigenvalues of nonsymmetric matrices can be arbitrarily large. If you have a matrix where left and right eigenvectors corresponding to the same eigenvalue are nearly orthogonal, a perturbation of size can move an eigenvalue by order /|l, r|, which can be enormous. There's no workaround for this other than recognizing that the eigenvalues are intrinsically sensitive and avoiding downstream operations that amplify the error.
Practical Code to Run Today
If you're using Python with NumPy, np.linalg.eigvals gives you eigenvalues directly. For a symmetric matrix, np.linalg.eigvalsh is faster and more accurate. In MATLAB, eigs is for sparse large matrices and eig is for dense ones. In C++, the Eigen library's SelfAdjointEigenSolver is optimized for symmetric problems and will beat the general solver by a wide margin on anything larger than 4x4. For a quick test, generate a random 5x5 symmetric matrix, compute its eigenvalues with both eigvalsh and eig, and compare. You'll see the symmetric solver converges faster and the results agree to machine precision while the general solver may show slightly more variance depending on the implementation. This is a normal difference, not a bug.
