近似: Padé 近似と有理補間

概要

有理近似は、対象関数を多項式の比 $R(x) = P(x)/Q(x)$ で近似する手法。 多項式近似 (Taylor 展開・Chebyshev 近似) が苦手とする場面で威力を発揮する:

• 極の表現: 分母 $Q(x)$ が零に近づくことで関数の発散・特異点を再現できる。多項式は有限値しか取れないため極を表せない。
• 漸近挙動: $x \to \infty$ で一定値に飽和する関数 ($\tanh$ など) を、分子・分母の最高次の比として自然に近似できる。
• 収束半径の超越: 元の Taylor 級数が発散する点でも、対角 Padé 近似はしばしば有効な値を返す (解析接続的なふるまい)。

同じ自由度 (係数の本数) を与えたとき、有理近似は多項式近似より広い範囲で誤差が小さくなることが多い。 代償として、近似分母の零点が偽の極 (および分子と対になった偽の零点 = Froissart doublet) を作る危険があり、極の位置を検査する規律が要る。

sangi の Approximation モジュールは、入力の種類に応じて 3 系統の有理近似を提供する:

• 単点の Taylor 係数から → padeApprox / padeTable (Padé 近似)
• Taylor 係数を連分数表現へ → taylorToSFraction / taylorToJFraction (qd アルゴリズム)
• 相異なる標本点の値から → rationalInterpolation (多点 Padé / Cauchy 補間)
• 行列指数 $e^A$ → matrixExpPade (スケーリング & スクェアリング)

関連 API: Approximation。 多項式・直交多項式系の近似は前ページ 最小二乗・Chebyshev・ミニマックス を参照。

Padé 近似

関数 $f$ の $x = 0$ まわりの Taylor 係数 $a_0, a_1, \ldots$ が与えられたとき、 $[M/N]$ Padé 近似は次の有理関数である:

$$R_{M/N}(x) = \frac{P_M(x)}{Q_N(x)} = \frac{p_0 + p_1 x + \cdots + p_M x^M}{1 + q_1 x + \cdots + q_N x^N}$$

分母は $Q_N(0) = 1$ と正規化する。係数は、$R_{M/N}$ の Taylor 展開が $f$ の Taylor 級数と 次数 $M+N$ まで一致するという条件で決まる:

$$f(x) - \frac{P_M(x)}{Q_N(x)} = O(x^{M+N+1})$$

係数マッチング条件と線形方程式系

両辺に $Q_N(x)$ を掛けて $x$ の各次数で係数を比較すると、次の条件が得られる (ここで $a_j = 0\ (j < 0)$ とする):

$$\sum_{l=0}^{N} q_l\, a_{m-l} = p_m \quad (0 \le m \le M), \qquad \sum_{l=0}^{N} q_l\, a_{m-l} = 0 \quad (M < m \le M+N)$$

$q_0 = 1$ なので、後半の $N$ 本は分母係数 $q_1, \ldots, q_N$ に関する 線形方程式系になる:

$$\begin{pmatrix} a_M & a_{M-1} & \cdots & a_{M-N+1} \\ a_{M+1} & a_M & \cdots & a_{M-N+2} \\ \vdots & & \ddots & \vdots \\ a_{M+N-1} & a_{M+N-2} & \cdots & a_M \end{pmatrix} \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \end{pmatrix} = -\begin{pmatrix} a_{M+1} \\ a_{M+2} \\ \vdots \\ a_{M+N} \end{pmatrix}$$

分母を解いたら、前半の関係から分子を畳み込みで陽に求める:

$$p_m = a_m + \sum_{l=1}^{\min(m,N)} q_l\, a_{m-l} \quad (0 \le m \le M)$$

sangi の padeApprox(taylor, M, N) はちょうどこの手順を踏む。 Toeplitz 構造の $N \times N$ 系を LU 分解 (algorithms::solve) で解き、 分子を畳み込みで構成して PadeResult<T> (numerator / denominator / valid) を返す。 分母系が特異 (退化 Padé) で解けなかったときは valid = false を返し、黙って誤った近似を返さない。 評価は evaluatePade が分子・分母をそれぞれ Horner 法で計算して比を取る。

Padé 表

$[M/N]$ を行 $M$・列 $N$ に並べた二次元配列をPadé 表と呼ぶ:

第 0 列 ($N=0$) は Taylor 多項式そのもの、対角 $[M/M]$ と 反対角 $M+N=\text{const}$ がそれぞれ収束の系列として重要。 実用では対角 Padé が最も安定して精度が高いことが多い。 sangi の padeTable(taylor, Mmax, Nmax) は $0 \le M \le M_{\max},\ 0 \le N \le N_{\max}$ の全セルを一括生成して vector<vector<PadeResult<T>>> として返す。 退化セルは valid = false のまま含める。各セルが $N \times N$ 系を解くため計算量は $O(M_{\max} \cdot N_{\max} \cdot N^2)$ で、対角だけが必要なら padeApprox を直接呼ぶ方が安い。

収束半径を超える近似

$f(x) = \log(1+x)$ や $f(x) = \tan x$ のように Taylor 級数の収束半径が有限な関数でも、 Padé 近似は収束円の外側で意味のある値を返すことがある。 分母の零点が $f$ の真の極の位置に近づくため、有理関数が解析接続の役割を果たすからである。 ただし、無関係な近接零点・極の対 (Froissart doublet) が現れることがあり、 高次の Padé では分母の零点を必ず点検する。

連分数と qd アルゴリズム

Padé 近似と表裏一体なのが連分数表現である。 Taylor 級数を連分数に変換しておくと、評価が桁落ちに強く、 途中で打ち切れば自動的に Padé 近似の系列が得られる。

Stieltjes 連分数 (S 分数)

S 分数は一段に一次の項を持つ形:

$$f(z) = \cfrac{c_0}{1 + \cfrac{a_1 z}{1 + \cfrac{a_2 z}{1 + \cfrac{a_3 z}{1 + \cdots}}}}$$

第 $k$ 近似 (途中で打ち切ったもの) が Padé 表のジグザグ系列に対応する。 sangi では StieltjesFraction<T> が先頭定数 b0 $= c_0$ と 係数列 a $= [a_1, a_2, \ldots]$、有限長で完結したかを示す terminated を保持する。

Jacobi 連分数 (J 分数)

J 分数は S 分数を二段ずつまとめた形で、二次の項を持つ:

$$f(z) = \cfrac{c_0}{1 - \beta_0 z - \cfrac{\alpha_1^2 z^2}{1 - \beta_1 z - \cfrac{\alpha_2^2 z^2}{1 - \beta_2 z - \cdots}}}$$

これは直交多項式の三項漸化式に対応し、第 $n$ 段までの J 分数は対角 $[n/n]$ Padé 近似に一致する。 sangi の JacobiFraction<T> は c0・係数列 beta ($\beta_0, \beta_1, \ldots$)・ alpha2 ($\alpha_1^2, \alpha_2^2, \ldots$) を保持する。

Taylor 級数 → 連分数変換 (qd アルゴリズム)

変換の核は Rutishauser の qd (quotient-difference) アルゴリズム。 Taylor 係数から $q$ 列と $e$ 列を交互に漸化式で生成する。第一行は

$$q_1^{(n)} = \frac{c_{n+1}}{c_n}, \qquad e_0^{(n)} = 0$$

から始め、菱形 (rhombus) 規則で表を埋めていく:

$$e_{k-1}^{(n)} = q_{k-1}^{(n+1)} - q_{k-1}^{(n)} + e_{k-2}^{(n+1)}, \qquad q_k^{(n)} = q_{k-1}^{(n+1)}\, \frac{e_{k-1}^{(n+1)}}{e_{k-1}^{(n)}}$$

得られた最上段の値から連分数の係数が決まる。S 分数は符号を反転して

$$a_{2k-1} = -q_k^{(0)}, \qquad a_{2k} = -e_k^{(0)},$$

J 分数は二段をまとめて

$$\beta_0 = q_1^{(0)}, \quad \beta_n = q_{n+1}^{(0)} + e_n^{(0)}\ (n \ge 1), \quad \alpha_n^2 = q_n^{(0)}\, e_n^{(0)}$$

で得る。sangi の taylorToSFraction(taylor, kmax) / taylorToJFraction(taylor, kmax) はこの qd 表を内部で構築する。 途中で $q_k$ または $e_k$ が $0$ になると (Padé の退化点)、そこで打ち切って terminated = true を立てる。先頭係数 $c_0 = 0$ の場合も $q_1 = c_1/c_0$ が定義できないため 即座に打ち切る。

連分数評価の安定性

連分数の評価は後退 (bottom-up) 法が安定: 最深段から始めて、上に向かって順に

$$\text{acc} \leftarrow 1 + \frac{a_i z}{\text{acc}}$$

を畳み上げ、最後に $f(z) = b_0 / \text{acc}$ とする (J 分数は二次項を含む同様の漸化)。 この方式は、分子・分母の多項式を別々に大きな値で評価してから比を取る素朴な方法に比べ、 中間値が $1$ 近傍に留まりやすく桁落ちに強い。sangi の evalSFraction / evalJFraction がこの後退評価を実装している。

有理補間 (Thiele / Cauchy)

Padé 近似は単一点の Taylor 係数に一致させる手法だった。 これに対し有理補間は、相異なる複数の標本点 $x_0, x_1, \ldots, x_{M+N}$ での値 $y_i = f(x_i)$ に一致する有理関数を求める。 多点 Padé・Newton-Padé・Cauchy 補間とも呼ばれる。

Thiele 型 (逆差分による連分数補間)

古典的な構成は Thiele の連分数補間で、Newton の差商の有理版である 逆差分 (reciprocal / inverse difference) $\varphi[\cdot]$ を用いて

$$R(x) = y_0 + \cfrac{x - x_0}{\varphi[x_0,x_1] + \cfrac{x - x_1}{\varphi[x_0,x_1,x_2] + \cfrac{x - x_2}{\ddots}}}$$

と連分数の形で表す。逆差分は

$$\varphi[x_i, x_{i+1}] = \frac{x_i - x_{i+1}}{y_i - y_{i+1}}, \qquad \varphi[x_0,\ldots,x_k] = \frac{x_{k-1} - x_k}{\varphi[x_0,\ldots,x_{k-2},x_{k-1}] - \varphi[x_0,\ldots,x_{k-2},x_k]}$$

の漸化で計算され、Padé における S 分数・J 分数の役割を多点で果たす。 Bulirsch-Stoer の常微分方程式ソルバや数値積分の外挿で、この有理外挿が使われる (標本を細かくして $h \to 0$ の極限を有理関数で外挿する)。

sangi の実装 — Cauchy 型線形系

sangi の rationalInterpolation(x, y, M, N) は、Thiele の逆差分を陽に組む代わりに、 補間条件を直接線形方程式系として解く Cauchy 型を採る。 各標本点で $P(x_i) - y_i\, Q(x_i) = 0$、すなわち

$$\sum_{k=0}^{M} p_k\, x_i^k \; - \; y_i \sum_{k=1}^{N} q_k\, x_i^k \; = \; y_i \qquad (i = 0, 1, \ldots, M+N)$$

を $q_0 = 1$ 正規化のもとで $M+N+1$ 元連立方程式として LU 分解で解く (標本点数はちょうど $M+N+1$ 個必要)。 係数列 numerator / denominator を RationalInterpResult<T> として返し、 評価は evalRational が Horner 法で行う。

Padé (単点) との違い・注意点

• 入力: Padé は 1 点の Taylor 係数、有理補間は相異なる多点の標本値。
• 標本の重複: $x_i$ に重複があると線形系が悪条件化する (本来は Hermite 型が必要)。 標本は相異なる点に取る。$x_i = 0$ に集約した Taylor 係数版が欲しいときは padeApprox を使う。
• 到達不能点: 一般に有理補間には、いかなる有理関数でも補間しきれない退化配置 (到達不能点) が存在する。 分母が標本点で零になる (極) 場合、分子も同時に零でない限り解は無効になる。

参考にした構成は Stoer & Bulirsch の有理補間 (§2.2.4) の枠組みである。

行列指数の Padé (スケーリング & スクェアリング)

Padé 近似はスカラー関数だけでなく行列値関数にも使える。 代表例が行列指数

$$e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!}$$

である。級数を直接打ち切る素朴な方法は、$\|A\|$ が大きいと項が一旦巨大化してから打ち消し合う (キャンセル) ため精度が出ない。実用的な標準解法がスケーリング & スクェアリング:

1. スケーリング: $\|A\|$ を見て $s$ を選び、$A/2^s$ のノルムを十分小さくする。
2. Padé: 小さくした $A/2^s$ に対角 Padé 近似 $e^{A/2^s} \approx R_{m/m}(A/2^s)$ を適用する。 指数関数の対角 Padé は分子・分母の係数が陽に与えられ、行列の解 (連立一次方程式) 1 回で評価できる。
3. スクェアリング: $e^A = \left(e^{A/2^s}\right)^{2^s}$ を $s$ 回の二乗で復元する。

sangi の matrixExpPade(A) は、ライブラリの行列指数ルーチン algorithms::expm への薄いラッパである。 expm はスケーリング & スクェアリングに $[13/13]$ 対角 Padé を中核として用い、 必要に応じてより高精度の Schur-Parlett 経路を持つ (Higham 2005 の方式)。引数 $m$ はソース互換のために残してあるが、 Padé の次数はルーチン内部で選ばれるため無視される。 正方行列でない入力は例外を投げる。

比較表

大まかな指針: Taylor 係数が手元にあって精度を引き上げたいなら Padé 近似、 その評価を安定化したい・打ち切りで Padé 系列が欲しいなら 連分数、 Taylor 係数でなく離散標本しかないなら 有理補間、 行列指数が必要なら スケーリング & スクェアリング。

参考文献

• Baker, G. A., & Graves-Morris, P. (1996). Padé Approximants. 2nd ed. Encyclopedia of Mathematics and Its Applications, Cambridge University Press.
• Cuyt, A., & Wuytack, L. (1987). Nonlinear Methods in Numerical Analysis. North-Holland.
• Rutishauser, H. (1957). Der Quotienten-Differenzen-Algorithmus. Birkhäuser.
• Stoer, J., & Bulirsch, R. (2002). Introduction to Numerical Analysis. 3rd ed. Springer. §2.2 (有理補間・Bulirsch-Stoer 外挿).
• Higham, N. J. (2005). "The scaling and squaring method for the matrix exponential revisited". SIAM Journal on Matrix Analysis and Applications, 26(4), 1179–1193.