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

技術要約：重心表示および連分数展開における導関数の安定な評価

問題提起
AAA法や貪欲なThiele連分数（TCF）選択アルゴリズムなどの登場により、有理関数近似への関心が再び高まっている。これらの手法は、$n$個のノードに基づく補間関数を生成し、$O(n)$の時間計算量で安定して評価できる。しかし、これらの表現形式における導関数の信頼性の高い計算には、大きな隔たりが存在する。

• 重心表示（Barycentric Representation）: SchneiderとWernerが$O(n)$の公式を提供しているものの、一般的な公式は補間ノード付近で減算による桁落ち（subtractive cancellation）が発生するため、数値的に不安定である。一方、ノードに対する代替公式は安定しているが、これとは異なる処理が必要となり、条件分岐を要する。
• TCF表示（TCF Representation）: 既存の1階導関数に関する$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）表示:

   • 安定性の証明: 著者らは、分子と分母のペアを用いるTCFの代替反復評価が、標準的な末尾順評価（tail-ordered evaluation）と同じ数値的安定性を持つことを証明する（定理1）。
   • 高階導関数: 前方モード自動微分を用いて、既存の1階導関数アルゴリズムを高階へと拡張する。著者らは、分子および分母多項式（$p_k$ および $q_k$）の導関数に関する漸化式（式21）を導出する。
   • アルゴリズム2: このアルゴリズムは、関数値とその $\mu$ 次までの導関数を $O(n)$ の時間で計算する。これには、中間値の再スケーリングによるアンダーフローへの実用的な対策が含まれている。

主要な結果
著者らはこれらのアルゴリズムを RationalFunctionApproximation.jl パッケージに実装し、区間 $[-1, 1]$ および単位円上の様々な関数（整関数、振動関数、メロモルフィック関数、および分岐点を持つ関数）を用いてテストを行った。

• 重心表示の安定性: 数値実験により、新しい公式が拡張精度（Double64）と一致する誤差を達成することが確認された。対照的に、元のSchneider–Wernerの汎用公式（式2）は、ノード付近で壊滅的な誤差（最大 $10^1$）を示すが、新しい公式は誤差をマシンイプシロン付近（$10^{-13}$ から $10^{-11}$）に維持する。
• TCFの安定性: TCF導関数アルゴリズムは、拡張精度における誤差とほぼ同一の誤差を倍精度で生成し、その実用的な安定性を裏付けている。
• 導関数の次数: これらの手法は、1階および2階導関数の計算に成功している。誤差の大きさは、一般に安定した重心表示の結果の1桁程度の範囲内に収まっている。
• 限界: 小さなパラメータ（$\epsilon \le 10^{-6}$）を持つ高度にパラメータ化された関数については、どちらの表現形式も一貫した正確なグローバル導関数を得られなかった。著者らは、これを導関数公式自体の不安定性ではなく、有理近似の微分における固有の不良条件設定（ill-conditioning）や、ドメイン付近に小さな残差を持つ極の存在に起因すると考えている。

意義と主張
本論文は、重心表示における導関数評価のための初の数値的に安定な手法（全階数の導関数に対する $O(n)$ アルゴリズムを含む）を提示していると主張している。TCF表示については、導関数アルゴリズムの高階への初の一般化を提供し、既存の1階導関数公式が実用的に安定していることを確認した。

著者らは控えめな主張を行っており、安定性の厳密な証明は欠いているものの、様々な領域における1階および2階導関数の数値実験によって、提案手法の堅牢性が確認されていると述べている。彼らは、有理近似の導関数の精度は、問題の条件付けによって根本的に制限されることが多く、関数近似自体の精度よりも低くなる傾向があることを強調している。本研究は、感度解析や最適化を必要とするアプリケーションにおける重要なボトルネックとなっている、安定した導関数評価のための実用的な解決策を提供するものである。