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

本文提出了一种利用迹、偏量不变量及判别式等四个不变量来数值稳定地计算实对角化 $3 \times 3$ 矩阵特征值的闭式方法，该方法不仅推导了误差界并验证了算法精度，还在保持精度的同时比 LAPACK 库快约十倍。

要点

这篇论文主要解决了一个数学和工程领域的“老难题”：如何快速、准确地算出 3x3 矩阵的特征值（Eigenvalues）。

为了让你轻松理解，我们可以把这篇论文的内容想象成一场**“寻找宝藏的探险”**，而那个“宝藏”就是矩阵的特征值。

1. 背景：为什么这是个难题？

想象你手里有一个神秘的 3x3 魔法盒子（矩阵），里面藏着三个数字（特征值）。在工程、物理和计算机图形学中，我们需要知道这三个数字来预测物体怎么变形、怎么旋转，或者光线怎么折射。

• 传统的做法（旧地图）： 以前，数学家们有一张“旧地图”（基于卡尔达诺和韦达的公式），告诉我们要怎么算出这三个数字。这张地图在大多数情况下很好用，就像在平坦的大路上开车。
• 遇到的问题（陷阱）： 但是，当这三个数字变得非常接近，甚至完全一样时（比如三个数字都是 1），旧地图就会失效。这时候，就像车子开进了沼泽地，稍微一点计算误差（就像车轮打滑），就会导致算出来的结果完全错误，甚至算出荒谬的数值。这在数学上叫“数值不稳定”。

2. 核心突破：绘制新地图（数值稳定的算法）

作者（Michal Habera 和 Andreas Zilian）决定重新绘制一张**“新地图”。他们不直接去算那个容易出错的“最终答案”，而是先计算四个关键的“路标”（不变量）**：

1. 迹 (Trace, $I_1$)： 就像盒子的“总重量”。
2. 偏量不变量 $J_2$ 和 $J_3$： 就像盒子的“形状特征”和“扭曲程度”。
3. 判别式 ($\Delta$)： 就像盒子的“稳定性指示器”，告诉我们这三个数字是分开得很远，还是挤在一起。

作者的魔法（新算法）：

• 旧方法的问题： 以前的算法在计算这些路标时，就像是用两个巨大的数字相减来得到一个小数字（例如 $1000000 - 999999 = 1$）。在计算机眼里，如果这两个大数字有一点点误差，结果就会变成$0$或者$2$，完全错了。这叫“灾难性抵消”。
• 新方法的高明之处： 作者设计了一套新的计算步骤（算法），就像换了一种更聪明的测量工具。他们利用**“对角线差值”**（比如 $A_{00} - A_{11}$）来代替直接的大数相减。
  • 比喻： 想象你要测量两栋楼的高度差。旧方法是分别量两栋楼（可能量不准），然后相减。新方法是直接拿一把尺子量两栋楼顶的高度差。这样，即使尺子有误差，对最终结果的影响也微乎其微。

3. 三大成果

A. 更精准（像激光一样准）

作者证明了他们的新算法在数学上是**“向后稳定”**的。

• 通俗解释： 这意味着，就算你的输入数据（那个魔法盒子）有一点点微小的抖动，新算法算出来的结果，依然非常接近那个“抖动后”的盒子的真实答案。
• 对比： 在测试中，当三个特征值几乎重合（最危险的情况）时，旧算法算出的结果可能错得离谱（误差大到 $10^{-8}$），而新算法的误差依然控制在计算机能达到的极限精度（$10^{-16}$），几乎完美。

B. 更快（像跑车一样快）

大家可能会想：“既然这么准，那肯定很慢吧？”

• 事实： 恰恰相反！新算法比目前工业界标准的“重型卡车”（LAPACK 库，一种非常成熟但复杂的迭代算法）快了10 倍。
• 比喻： LAPACK 就像一辆装备精良的坦克，什么地形都能走，但启动慢、转弯慢。作者的新算法就像一辆F1 赛车，在特定的赛道（3x3 矩阵）上，它既快又稳。

C. 更实用（拯救了工程师）

论文最后展示了一个实际案例：莫尔 - 库仑屈服函数（用于预测土壤或岩石何时会崩塌）。

• 场景： 在计算土壤会不会塌方时，如果算错了特征值，工程师可能会误判，导致桥梁或大坝设计不安全。
• 结果： 使用旧算法，在土壤处于临界状态（特征值重合）时，计算出的“安全系数”会乱跳，导致误报。使用作者的新算法，无论土壤状态多复杂，计算结果都稳稳当当，始终保持在机器精度范围内。

4. 总结：这篇论文意味着什么？

这就好比作者发现了一种**“防抖动的数学滤镜”**。

• 以前： 算 3x3 矩阵的特征值，就像在暴风雨中走钢丝，一旦遇到特殊情况（数字重合），很容易掉下去。
• 现在： 作者给这根钢丝加了护栏（新算法），并且铺上了防滑垫。
  • 对于大多数情况（特征值分开），它比现有的标准方法快 10 倍。
  • 对于最危险的情况（特征值重合），它依然能保持极高的精度，而旧方法会彻底失效。

虽然对于极度混乱的“沼泽地”（极度病态的矩阵），传统的“重型坦克”（LAPACK）可能还是更可靠，但对于绝大多数工程应用，作者提供的这套**“新地图”**是更快、更准、更安全的最佳选择。

一句话总结： 作者发明了一套新的数学公式，让计算机在计算 3x3 矩阵的“核心数字”时，既快如闪电，又稳如泰山，彻底解决了老方法在特殊情况下容易“翻车”的难题。

技术摘要

这是一份关于3×3 矩阵特征值闭式表达式的数值稳定评估的详细技术总结。该论文由 Michal Habera 和 Andreas Zilian 撰写，旨在解决传统三角函数公式在特征值重根（coalescing eigenvalues）情况下数值不稳定的问题。

1. 研究背景与问题 (Problem)

• 核心问题：对于具有实谱的可对角化 $3 \times 3$矩阵，传统的特征值闭式解（基于 Cardano 和 Viète 的三次方程求根公式，利用矩阵迹$I_1$和偏量不变量$J_2, J_3$ 计算）在有限精度算术下存在严重的数值不稳定性。
• 具体痛点：
  • 当特征值重合（重根）或接近重合时，传统的 $\arccos$ 公式在 $x \to 1$ 附近会导致巨大的舍入误差。
  • 判别式 $\Delta = 4J_2^3 - 27J_3^2$ 的计算涉及两个相近大数的相减，极易发生灾难性抵消（catastrophic cancellation）。
  • 现有的文献往往缺乏严格的误差分析，且许多“数值稳定”的算法仅针对对称矩阵，无法推广到非对称矩阵。
• 现有替代方案的局限：虽然迭代算法（如 LAPACK 中的 QR 算法）通常稳定，但它们计算量大，且不支持符号微分（在优化和机器学习应用中至关重要），也无法像闭式解那样直接保证特征值的排序。

2. 方法论 (Methodology)

作者提出了一套基于四个关键不变量的数值稳定闭式评估算法，并进行了严格的条件数分析和前向误差界推导：

1. 不变量定义：

   • $I_1$：矩阵的迹 (Trace)。
   • $J_2, J_3$：偏量不变量 (Deviatoric invariants)。
   • $\Delta$：判别式 (Discriminant)。

2. 核心算法改进：

   • $I_1$：直接求和，证明其具有前向稳定性。
   • $J_2$：提出了一种基于对角线差值（diagonal differences）和非对角线乘积的算法（Algorithm 2）。通过先计算 $A_{ii} - A_{jj}$，避免了直接计算偏量矩阵时的舍入误差，确保对于缩放单位矩阵 $A = \alpha I$，计算结果精确为零。
   • $J_3$：提出了基于对角线差值、非对角线乘积和混合项的展开公式（Algorithm 5），推广了 Harari 和 Albocher 针对对称矩阵的工作至非对称矩阵。
   • $\Delta$：利用 Parlett 的因子分解思想，采用乘积和公式（Sum-of-products formula，Algorithm 8）。该公式将判别式表示为 14 项的加权和，其中每一项都包含对角线差值。当矩阵接近重根时，这些因子趋于零，从而避免了直接计算 $4J_2^3 - 27J_3^2$ 时的抵消问题。
   • 特征值计算：利用改进后的不变量，结合 $\arctan$ 公式（替代 $\arccos$）计算三重角 $\phi$，进而通过三角公式求解特征值（Algorithm 9）。

3. 理论分析：

   • 定义了相对偏量后向稳定性（Relative deviatoric backward stability），即扰动 $\delta A$ 的偏量部分 $\|\text{dev}(\delta A)\|$ 相对于 $\|\text{dev}(A)\|$ 是小的。
   • 利用 Bauer-Fike 定理推导了特征值计算的前向误差界，指出误差与特征基的条件数 $\kappa_2(U)$ 成正比。

3. 关键贡献 (Key Contributions)

1. 严格的误差分析：填补了现有文献中缺乏对闭式特征值算法进行严格前向/后向误差分析的空白。推导了 $I_1, J_2, J_3, \Delta$ 的误差上界。
2. 通用性：提出的算法不仅适用于对称矩阵，还适用于非对称且可对角化的 $3 \times 3$ 矩阵。
3. 针对 $J_2 \to 0$ 的优化：特别解决了当 $J_2$ 趋近于 0（即矩阵接近单位矩阵的倍数）时的数值稳定性问题，证明了新算法在此极限情况下的准确性。
4. 开源库：发布了开源库 eig3x3，包含 C11 实现（带 Python 接口），提供了从朴素算法到稳定算法的完整对比。

4. 实验结果 (Results)

作者通过基准测试（Benchmark）对比了三种方法：朴素多项式求和（Naive）、基于 NumPy 张量运算的朴素方法（Naive Tensor）和本文提出的稳定算法（Present），并与 LAPACK 的 DGEEV 进行了对比。

• 数值稳定性：
  • $J_2$ 和 $J_3$：本文算法在所有测试案例（包括特征基条件数良好和病态的情况）中均满足理论推导的稳定性界。特别是在 $J_2 \to 0$ 的极限情况下，误差随 $\delta^2$ 减小，达到了机器精度级别。
  • 判别式 $\Delta$：对于特征基条件数良好（well-conditioned）或正交的矩阵，算法表现优异；但在特征基极度病态（ill-conditioned）的情况下，精度会有所下降，但仍优于朴素方法。
  • 特征值：对于特征基条件数良好的矩阵，本文算法的误差与 LAPACK 相当，且满足稳定性界。对于病态矩阵，LAPACK 表现略好，但本文算法仍显著优于朴素闭式解。
• 性能：
  • 在具有接近重根的挑战性测试用例中，本文提出的闭式算法比高度优化的 LAPACK 库（DGEEV）快约 10 倍（平均耗时 38ns vs 396ns）。
  • 由于是闭式解，代码可内联（inline），无函数调用开销，在性能敏感应用中优势更大。
• 实际应用验证：
  • 在Mohr-Coulomb 屈服函数（岩土工程中预测材料失效）的计算中，朴素算法在特征值重合附近（Lode 角特定值）产生高达 $10^{-8}$的误差，而本文算法在整个应力状态范围内均保持接近机器精度（$10^{-16}$）的误差。

5. 意义与结论 (Significance & Conclusion)

• 理论价值：证明了闭式解在特定条件下（特征基条件数良好）可以达到与迭代法相当的数值稳定性，同时提供了严格的误差界限。
• 工程价值：
  • 高性能：比标准库快一个数量级，适用于实时计算、嵌入式系统或高频计算场景。
  • 符号微分友好：闭式解允许直接进行符号微分，这对于优化问题、机器学习中的梯度计算以及计算塑性力学中的返回映射算法（return mapping algorithms）至关重要。
  • 特征值排序：闭式解天然保证特征值有序（$\lambda_1 \le \lambda_2 \le \lambda_3$），这在主应力/主应变分析中非常重要，而迭代法通常不保证排序。
• 局限性：对于特征基极度病态（ill-conditioned eigenbasis）的矩阵，迭代法（如 LAPACK）仍然是精度最高的选择。未来的工作将致力于进一步证明 $J_3$ 和 $\Delta$ 在病态情况下的前向稳定性，并进行针对特定 CPU 架构的 SIMD 优化。

总结：该论文提供了一套经过严格数学证明、数值稳定且高效的 $3 \times 3$ 矩阵特征值闭式解算法，解决了传统方法在重根附近的数值崩溃问题，并在保持高精度的同时显著提升了计算速度，具有广泛的工程应用前景。