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

Technical Summary: Stable Evaluation of Derivatives for Barycentric and Continued Fraction Representations

Problem Statement
Rational function approximation has seen renewed interest due to algorithms like AAA and greedy Thiele continued fraction (TCF) selection, which produce interpolants on $n$ nodes that can be evaluated stably in $O(n)$ time. However, a significant gap exists in the reliable computation of derivatives for these representations.

• Barycentric Representation: While Schneider and Werner provided $O(n)$ formulas for derivatives, the generic formula is numerically unstable near interpolation nodes due to subtractive cancellation. The alternative formula for nodes is stable but distinct, requiring conditional logic.
• TCF Representation: An existing $O(n)$ iterative formula for the first derivative (based on forward-mode automatic differentiation) lacked a formal investigation into its numerical stability, and no generalization to higher-order derivatives had been established.

Methodology
The authors propose new algorithms to evaluate derivatives for both representations while maintaining numerical stability and computational efficiency.

1. Barycentric Representation:

   • Analysis of Instability: The authors demonstrate that the standard formula $r'(z) = [N'(z) - r(z)D'(z)]/D(z)$ suffers from dominant subtractive cancellation when $z$ is close to a node $z_j$.
   • Double-Sum Solution: They derive a double-summation formula (Equation 8) that eliminates this cancellation. While stable, this approach has $O(n^2)$ complexity.
   • Stable $O(n)$ Algorithm: To retain linear complexity, the authors introduce a modified approach. By defining $\tilde{N}(z) = \epsilon N(z)$ and $\tilde{D}(z) = \epsilon D(z)$ where $\epsilon = z - z_j$ (with $z_j$ being the nearest node), they derive a stable recurrence relation (Equation 13) for derivatives of all orders. This method avoids the cancellation by working with the scaled numerator and denominator.
   • Algorithm 1: The proposed method identifies the nearest node in $O(n)$ time and applies the modified recurrence to compute derivatives up to order $\mu$.

2. Thiele Continued Fraction (TCF) Representation:

   • Stability Proof: The authors prove (Theorem 1) that the alternative iteration for evaluating TCFs (using numerator/denominator pairs) has the same numerical stability as the standard tail-ordered evaluation.
   • Higher-Order Derivatives: They extend the existing first-derivative algorithm to arbitrary orders using forward-mode automatic differentiation. They derive recurrence relations (Equation 21) for the derivatives of the numerator and denominator polynomials ($p_k$ and $q_k$).
   • Algorithm 2: This algorithm computes the function value and its derivatives up to order $\mu$ in $O(n)$ time. It includes a practical safeguard against numerical underflow by rescaling intermediate values.

Key Results
The authors implemented these algorithms in the RationalFunctionApproximation.jl package and tested them on various functions (entire, oscillatory, meromorphic, and those with branch points) over the interval $[-1, 1]$ and the unit circle.

• Barycentric Stability: Numerical experiments confirm that the new formula achieves errors consistent with extended-precision arithmetic (Double64). In contrast, the original Schneider–Werner generic formula (Equation 2) exhibits catastrophic errors (up to $10^1$) near nodes, while the new formula maintains errors near machine epsilon ($10^{-13}$ to $10^{-11}$).
• TCF Stability: The TCF derivative algorithm produces errors in double precision that are nearly identical to those in extended precision, confirming its practical stability.
• Order of Derivatives: The methods successfully compute first and second derivatives. The error magnitudes are generally within an order of magnitude of the stable barycentric results.
• Limitations: For highly parameterized functions with small parameters ($\epsilon \le 10^{-6}$), neither representation consistently yields accurate global derivatives. The authors attribute this to the inherent ill-conditioning of differentiating rational approximations and the potential presence of poles with small residuals near the domain, rather than instability in the derivative formulas themselves.

Significance and Claims
The paper claims to present the first numerically stable methods for evaluating derivatives in the barycentric representation, including an $O(n)$ algorithm for derivatives of all orders. For the TCF representation, the work provides the first generalization of the derivative algorithm to higher orders and confirms that the existing first-derivative formula is stable in practice.

The authors are modest in their claims, noting that while they lack formal proofs of stability, the numerical experiments on first and second derivatives across various domains confirm the robustness of the proposed methods. They emphasize that the accuracy of rational derivative approximations is fundamentally limited by the conditioning of the problem, often yielding lower accuracy than the function approximation itself. The work provides a practical solution for stable derivative evaluation in rational approximation workflows, addressing a critical bottleneck in applications requiring sensitivity analysis or optimization.