Numerically stable evaluation of closed-form expressions for eigenvalues of $3 \times 3$ matrices

This paper presents numerically stable, closed-form algorithms for computing the eigenvalues of real, diagonalizable $3 \times 3$ matrices using four specific invariants, demonstrating that these methods are both more accurate for repeated eigenvalues and approximately ten times faster than the LAPACK library while maintaining comparable precision.

The Gist

Imagine you have a magical box (a 3x3 matrix) that holds three hidden numbers inside it, called eigenvalues. These numbers are crucial for engineers and scientists because they tell us how a material will break, how a fluid will swirl, or how a bridge will vibrate.

For decades, mathematicians have had a "recipe" (a closed-form formula) to find these hidden numbers instantly, without needing to guess and check. However, this recipe has a fatal flaw: it's like a house of cards. If the hidden numbers are very close to each other (or identical), the recipe collapses, and the computer spits out garbage results.

This paper introduces a reinforced, earthquake-proof version of that recipe. Here is the breakdown of how they did it, using simple analogies.

1. The Problem: The "Catastrophic Cancellation" Trap

Think of the old recipe as trying to measure the height of a mountain by subtracting two huge numbers that are almost identical.

• Old Way: You have a number like 1,000,000.000001 and you subtract 1,000,000.000000.
• The Result: Mathematically, the answer is 0.000001. But in a computer's brain (which has limited precision), those tiny decimals get lost in the noise. The computer might think the answer is 0 or 0.000002.
• The Consequence: When the hidden numbers in the matrix are close together (like a mountain range with peaks of similar height), the old formula gets confused and produces wildly inaccurate results.

2. The Solution: Four "Stable" Ingredients

The authors realized that instead of trying to solve the whole puzzle at once, they should break it down into four specific ingredients (called invariants) and calculate each one with extreme care.

Think of these ingredients as the foundation, the shape, the weight, and the balance of a building:

1. Trace ($I_1$): The total "weight" of the matrix. (Easy to calculate, very stable).
2. Deviatoric Invariant 2 ($J_2$): A measure of how "stretched" or "squashed" the matrix is.
3. Deviatoric Invariant 3 ($J_3$): A measure of the "twist" or "handedness" of the matrix.
4. Discriminant ($\Delta$): A special number that tells us if the hidden numbers are about to merge (coalesce).

The Innovation:
The authors didn't just use the standard formulas for these ingredients. They rewrote the math for $J_2$, $J_3$, and $\Delta$ to avoid the "subtraction trap."

• Analogy: Instead of measuring the difference between two tall buildings directly (which is prone to error), they measured the difference between the roofs and the foundations separately and combined the results in a way that cancels out the noise. They used a "sum-of-products" method, which is like building a wall with interlocking bricks rather than stacking loose stones.

3. The "Well-Conditioned" vs. "Ill-Conditioned" Test

The paper tests their new recipe against two types of scenarios:

• The "Well-Behaved" Case: The hidden numbers are distinct, or the matrix is perfectly symmetrical (like a perfect sphere).
  • Result: The new recipe is perfect. It hits the target every time.
• The "Chaos" Case: The hidden numbers are merging, and the matrix is twisted (ill-conditioned).
  • Result: The new recipe is still very good, but it hits a limit. In these extreme chaos cases, the "gold standard" iterative methods (like the LAPACK library used by supercomputers) are still slightly more accurate. However, the new recipe is much closer to the truth than the old broken recipe.

4. Speed: The Sports Car vs. The Tank

Why not just use the "gold standard" (LAPACK) for everything?

• LAPACK is like a heavy, armored tank. It is incredibly accurate and can handle any terrain, but it is slow and heavy. It takes many steps to get to the answer.
• The New Recipe is like a Formula 1 sports car. It is a direct, closed-form path. It doesn't need to take detours.
• The Result: The new recipe is 10 times faster than the tank. In applications where you need to calculate these numbers millions of times per second (like simulating a car crash or training an AI), that speed difference is massive.

5. Real-World Impact: The "Mohr-Coulomb" Test

To prove it works in the real world, the authors tested it on a model used by civil engineers to predict when soil or rock will crumble under pressure (the Mohr-Coulomb yield function).

• Old Recipe: Near the point where the soil is about to fail, the old recipe got the math wrong, potentially leading engineers to think a dam is safe when it's actually about to collapse.
• New Recipe: It calculated the failure point with near-perfect precision, even when the stress levels were critical.

Summary

This paper is about fixing a broken calculator.

• The Problem: The old way of calculating eigenvalues crashes when numbers get too similar.
• The Fix: A new, mathematically robust way to calculate the intermediate steps, avoiding the "noise" of computer arithmetic.
• The Benefit: It is 10x faster than the current industry standard and much more accurate in critical situations, making it a game-changer for engineering simulations, physics engines, and machine learning.

It's the difference between using a shaky ruler that breaks when you measure two identical lines, versus a laser measure that stays precise no matter how close the lines get.

Technical Summary

Here is a detailed technical summary of the paper "Numerically Stable Evaluation of Closed-Form Expressions for Eigenvalues of 3 × 3 Matrices" by Michal Habera and Andreas Zilian.

1. Problem Statement

The computation of eigenvalues for $3 \times 3$ real matrices is a fundamental task in engineering, physics, and computer science. While iterative methods (like the QR algorithm in LAPACK) are robust, closed-form expressions are often preferred for performance-critical applications, symbolic differentiation, and scenarios requiring ordered eigenvalues.

However, classical closed-form solutions based on Cardano's and Viète's trigonometric formulas suffer from numerical instability, particularly when:

• Eigenvalues are repeated or coalesce (e.g., $J_2 \to 0$ or $\Delta \to 0$).
• The matrix has an ill-conditioned eigenbasis (non-normal matrices).
• Standard formulas involve catastrophic cancellation (subtracting nearly equal large numbers) when computing invariants like the discriminant $\Delta = 4J_2^3 - 27J_3^2$.

Existing literature often lacks rigorous forward error analysis, and many "stable" algorithms fail to generalize to non-symmetric matrices or degrade significantly when the eigenbasis is ill-conditioned.

2. Methodology

The authors propose a framework for evaluating eigenvalues using four key matrix invariants:

1. Trace ($I_1$): $tr(A)$.
2. Deviatoric Invariant ($J_2$): Related to the second invariant of the deviatoric part.
3. Deviatoric Invariant ($J_3$): The determinant of the deviatoric part.
4. Discriminant ($\Delta$): $4J_2^3 - 27J_3^2$.

The methodology involves three main pillars:

A. Rigorous Error Analysis

The paper adopts the standard IEEE 754 floating-point model and defines backward stability and forward error bounds. The authors derive tight theoretical bounds for the absolute forward error of each invariant, relating them to the condition number of the matrix eigenbasis ($\kappa_2(U)$).

B. Development of Stable Algorithms

Instead of using naive monomial expansions or standard tensor operations (which accumulate rounding errors), the authors propose specialized algorithms:

• For $I_1$: Simple summation (trivially stable).
• For $J_2$: An algorithm (Algorithm 2) that computes $J_2$ using diagonal differences ($A_{ii} - A_{jj}$) and off-diagonal products. This ensures that for a scaled identity matrix ($A = \alpha I$), the result is exactly zero, avoiding catastrophic cancellation near $J_2 = 0$.
• For $J_3$: An expansion (Algorithm 5) based on diagonal differences and mixed products, generalizing previous work for symmetric matrices to non-symmetric cases.
• For $\Delta$: A sum-of-products formula (Algorithm 8) derived from the Cauchy–Binet formula. This expresses the discriminant as a sum of 14 terms where factors vanish as the matrix approaches multiple eigenvalues, preventing the subtraction of large, nearly equal numbers inherent in the standard $4J_2^3 - 27J_3^2$ formula.

C. Eigenvalue Reconstruction

The eigenvalues are reconstructed using the trigonometric formula:
$$\lambda_k = \frac{1}{3} \left( I_1 + 2\sqrt{3J_2} \cos\left(\frac{\phi + 2\pi k}{3}\right) \right)$$
Crucially, the triple-angle $\phi$ is computed using arctan (Eq. 4) rather than arccos to improve stability near repeated eigenvalues.

3. Key Contributions

1. Rigorous Forward Error Bounds: The paper provides the first rigorous derivation of forward error bounds for $J_2$, $J_3$, and $\Delta$ in the context of non-symmetric matrices. It introduces the concept of relative deviatoric backward stability, showing that controlling the error in the deviatoric part of the perturbation is sufficient for accuracy.
2. Stable Algorithms for Non-Symmetric Matrices: Unlike previous works limited to symmetric matrices, the proposed algorithms (Algorithms 2, 5, and 8) are designed for general real, diagonalizable $3 \times 3$ matrices.
3. Performance vs. Accuracy Trade-off Analysis: The authors demonstrate that while iterative solvers (LAPACK) are universally stable, the proposed closed-form method is ~10x faster while maintaining machine precision accuracy for matrices with well-conditioned eigenbases.
4. Open-Source Implementation: The authors released the eig3x3 library (C11 with Python wrappers), providing production-ready implementations of these stable algorithms.

4. Results

The authors benchmarked their algorithms against "naive" implementations and the highly optimized LAPACK routine DGEEV using two specific test paths:

• Path 1 ($D_1$): Approaching a triple eigenvalue ($J_2 \to 0, J_3 \to 0$).
• Path 2 ($D_2$): Approaching a double eigenvalue ($\Delta \to 0$).

Key Findings:

• $J_2$ Algorithm: The proposed algorithm is accurate and stable for all cases, including those with ill-conditioned eigenbases ($\kappa_2 \approx 9000$). It satisfies the theoretical error bounds.
• $J_3$ and $\Delta$ Algorithms: These are stable for well-conditioned and orthogonal eigenbases. However, for ill-conditioned eigenbases, the error degrades and exceeds the theoretical bounds, though it remains significantly better than naive methods.
• Eigenvalue Computation:
  • Naive methods: Produced errors up to $10^{-5}$(for triple eigenvalues) and$10^{-8}$ (for double eigenvalues).
  • Proposed Method: Achieved errors close to machine precision ($\approx 10^{-16}$) for well-conditioned matrices.
  • LAPACK DGEEV: Consistently achieved machine precision across all cases, including ill-conditioned ones, but was significantly slower.
• Performance: The proposed algorithm was approximately 10 times faster than LAPACK DGEEV (38 ns vs. 396 ns per evaluation) on a challenging test case.

5. Significance and Applications

The paper bridges the gap between theoretical numerical stability and practical engineering needs:

• Computational Plasticity: The authors demonstrated the utility of their method by evaluating the Mohr–Coulomb yield function. Naive methods produced massive errors ($10^{-8}$) near yield surfaces where eigenvalues coalesce, potentially leading to incorrect failure predictions. The stable algorithm maintained machine precision accuracy.
• Symbolic Differentiation: Because the solution is closed-form, it allows for exact symbolic differentiation (e.g., computing eigenprojectors), which is essential for optimization and machine learning applications involving matrix functions.
• Real-Time Systems: The 10x speedup makes these algorithms viable for real-time simulations (e.g., fluid dynamics, structural analysis) where iterative solvers are too slow, provided the matrices are not extremely ill-conditioned.

Conclusion

The authors successfully developed a numerically stable, closed-form evaluation method for $3 \times 3$ matrix eigenvalues. While iterative solvers remain the gold standard for the most ill-conditioned problems, this work provides a robust, high-performance alternative for the vast majority of practical engineering cases, backed by rigorous error analysis and open-source software.