Stable evaluation of derivatives for barycentric and continued fraction representations of rational functions

技术摘要：重心表示与 Thiele 连分数表示的稳定导数求值

问题陈述
由于 AAA 算法和贪婪 Thiele 连分数（TCF）选择算法等技术的出现，有理函数逼近重新受到了关注，这些算法生成的插值在 $n$ 个节点上可以实现 $O(n)$ 时间内的稳定求值。然而，在这些表示法中可靠地计算导数仍存在显著空白。

• 重心表示（Barycentric Representation）： 虽然 Schneider 和 Werner 提供了 $O(n)$ 阶导数的公式，但通用公式在靠近插值节点时会因减法抵消（subtractive cancellation）而导致数值不稳定。另一种针对节点的替代公式虽然稳定，但与通用公式不同，需要条件逻辑。
• TCF 表示： 现有的基于前向模式自动微分的 $O(n)$ 一阶导数迭代公式缺乏对其数值稳定性的正式研究，且尚未建立向高阶导数推广的方法。

方法论
作者提出了新的算法，用于在保持数值稳定性和计算效率的同时，为这两种表示法求取导数。

1. 重心表示：

   • 不稳定性分析： 作者证明了标准公式 $r'(z) = [N'(z) - r(z)D'(z)]/D(z)$ 在 $z$ 靠近节点 $z_j$ 时会遭受严重的减法抵消。
   • 双重求和解法： 他们推导出了一个消除这种抵消的双重求和公式（公式 8）。虽然该方法是稳定的，但其复杂度为 $O(n^2)$。
   • 稳定的 $O(n)$ 算法： 为了保留线性复杂度，作者引入了一种改进方法。通过定义 $\tilde{N}(z) = \epsilon N(z)$ 和 $\tilde{D}(z) = \epsilon D(z)$，其中 $\epsilon = z - z_j$（$z_j$ 为最近的节点），他们推导出了一个用于计算各阶导数的稳定递推关系（公式 13）。该方法通过处理缩放后的分子和分母来避免抵消。
   • 算法 1： 该方法在 $O(n)$ 时间内识别最近的节点，并应用改进后的递推关系来计算最高 $\mu$ 阶的导数。

2. Thiele 连分数（TCF）表示：

   • 稳定性证明： 作者证明（定理 1），使用分子/分母对进行评估的 TCF 替代迭代法具有与标准尾部排序评估相同的数值稳定性。
   • 高阶导数： 他们利用前向模式自动微分，将现有的一阶导数算法扩展到任意阶。他们推导了用于计算分子和分母多项式（$p_k$ 和 $q_k$）导数的递推关系（公式 21）。
   • 算法 2： 该算法在 $O(n)$ 时间内计算函数值及其最高 $\mu$ 阶导数。它包含一个用于防止数值下溢的实用保护机制，通过对中间值进行重缩放来实现。

关键结果
作者在 RationalFunctionApproximation.jl 软件包中实现了这些算法，并在 $[-1, 1]$ 区间及单位圆上针对各种函数（整函数、振荡函数、亚纯函数以及带有分歧点的函数）进行了测试。

• 重心表示稳定性： 数值实验证实，新公式实现的误差与扩展精度算术（Double64）一致。相比之下，原始的 Schneider–Werner 通用公式（公式 2）在靠近节点时会出现灾难性误差（高达 $10^1$），而新公式能将误差维持在接近机器 epsilon 的水平（$10^{-13}$ 到 $10^{-11}$）。
• TCF 稳定性： TCF 导数算法在双精度下产生的误差与扩展精度下的误差几乎相同，证实了其在实际中的稳定性。
• 导数阶数： 这些方法成功计算了一阶和二阶导数。其误差量级通常在稳定重心表示结果的一个数量级之内。
• 局限性： 对于具有微小参数（$\epsilon \le 10^{-6}$）的高度参数化函数，两种表示法都无法一致地获得准确的全局导数。作者将其归因于有理逼近求导过程中固有的病态性，以及可能存在的残差极小的极点位于定义域附近，而非导数公式本身的稳定性问题。

意义与主张
本文声称提出了首个在重心表示中实现数值稳定求导的方法，包括一个用于计算各阶导数的 $O(n)$ 算法。对于 TCF 表示，这项工作提供了首次向高阶导数的推广，并确认了现有的一阶导数公式在实践中是稳定的。

作者的声明较为谦逊，指出虽然他们缺乏形式化的稳定性证明，但在各种定义域上对一阶和二阶导数的数值实验证实了所提方法的鲁棒性。他们强调，有理导数逼近的精度从根本上受限于问题的条件数，通常比函数逼近本身的精度要低。这项工作为需要进行灵敏度分析或优化的应用场景提供了一个稳定的导数计算实用解决方案，解决了有理逼近工作流中的一个关键瓶颈。