Eig3x3.Eigenvalues — the Habera–Zilian invariant pipeline

Given a real symmetric matrix A to decompose, performs Habera-Zilian's method [HZ25] to compute an ordered vector of eigenvalues, Eigval3.

Algorithm

The eigenvalues of A are the roots of its cubic characteristic polynomial. While cubics have a closed-form trigonometric solution, evaluating it naively loses accuracy exactly when two eigenvalues are close.

The solution is to sidestep approaching the polynomial directly, and instead recenter by the mean of the eigenvalues, which the trace i1 provides trivially, so the remaining unknowns sum to zero. Then, measure the spread of the recentered eigenvalues with two moments — j2, roughly their variance, and j3, roughly their skewness — computed from differences of diagonal entries so that nearly-equal eigenvalues do not cancel away.

The discriminant, Δ = 4J₂³ − 27J₃², (delta) is evaluated as a sum of squares so no subtraction of nearly-equal eigenvalues can occur and cause floating-point instability. While the classical formula uses the arccos, the required angle can be expressed as the arctan2, leveraging its numerical stability near zero. Finally, the eigenvalues are computed λₖ = (I₁ + 2√(3J₂)·cos(φ/3 + 2πk/3))/3 for k = { 1, 2, 3 } (eigvals).

Provenance

Habera & Zilian, "Numerically stable evaluation of closed-form expressions for eigenvalues of 3×3 matrices", arXiv:2511.00292v2 (2025). [HZ25] (CC BY 4.0).

Specifically:

• invariants I₁ (Alg. 1), J₂ (Alg. 2), J₃ (Alg. 5),
• the Algorithm 8 sum-of-squares discriminant,
• verified Algorithm 8 through factorization of the original Habera-Zilian algorithm ([HZ21], Eq. 29),
• quadrant-safe angle φ = atan2(√(27Δ), 27J₃) (Eq. 4),
• and ordered eigenvalues λ₁ ≤ λ₂ ≤ λ₃ (Eq. 2).

Deviations:

• Eq. 4's arctan of the ratio is realized as atan2 of numerator and denominator, which is identical when J₃ > 0, quadrant-correct when J₃ < 0, and NaN-free at J₃ = Δ = 0.

Visibility

Exposes eigvals.

The rest of this module is internal, package-private.

def Eig3x3.eigvals (A : SymmMat3) : Eigval3

Eigenvalues in increasing order ([HZ25] Eq. 2 with the Eq. 4 arctan angle). atan2 yields φ ∈ [0, π]; k = 1, 2, 3 then gives λ₁ ≤ λ₂ ≤ λ₃. For a scaled identity, J₂ = J₃ = Δ = 0, φ = atan2(0,0) = 0, and all three eigenvalues come out as exactly I₁/3.

Postcondition (contract): the result satisfies l₀ ≤ l₁ ≤ l₂ enforced by the final sort.

Equations

• One or more equations did not get rendered due to their size.