補間: 区分・形状保存補間

概要

補間とは、与えられた節点 $(x_i, y_i)$ ($i = 0, \ldots, n$) をちょうど通る関数を構成し、節点の間の値を推定する操作。 最も素朴な方法は、$n+1$ 点をすべて通る次数 $n$ の多項式補間を 1 本で構成することだが、これは次数が上がると破綻する。

Runge 現象

等間隔の節点で高次の多項式補間を行うと、区間の端付近で補間多項式が激しく振動し、節点を増やすほど振動が悪化する。 古典的な例は $f(x) = 1/(1 + 25 x^2)$ を区間 $[-1, 1]$ で等間隔補間したもので、端で誤差が指数的に発散する。これを Runge 現象と呼ぶ。

原因は、等間隔節点におけるラグランジュ基底 $\ell_i(x)$ の振幅が端で爆発的に大きくなることにある。 対策は大きく 2 つ:

• 節点を端に密にする (Chebyshev 節点など) — 1 本の多項式のまま振動を抑える。
• 区分・局所的にする — 各区間を低次 (1〜3 次) の多項式でつなぎ、誤差を局所化する。これが区分補間。

区分補間は、評価点が属する 1 区間の係数だけで値が決まる局所性を持つ。 節点を 1 つ動かしても影響が近傍に限られるため、データ駆動の応用で扱いやすい。 本ページでは、sangi が実装する区分・形状保存系の手法 — 3 次スプライン、PCHIP、修正 Akima、Catmull-Rom、5 次 Hermite — を、その基礎となるラグランジュ・Newton 補間とあわせて解説する。

関連 API: Interpolation。

ラグランジュ補間と Newton 分割差分

区分手法に入る前に、その土台となる大域多項式補間の 2 つの古典形式を整理する。 どちらも同じ補間多項式 $p(x)$ を表すが、構成と評価の効率が異なる。

ラグランジュ基底

節点 $x_0, \ldots, x_n$ に対し、ラグランジュ基底を

$$\ell_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^{n} \frac{x - x_j}{x_i - x_j}$$

と定めると、$\ell_i(x_k) = \delta_{ik}$ ($i = k$ で 1、それ以外で 0) が成り立つ。補間多項式は

$$p(x) = \sum_{i=0}^{n} y_i\, \ell_i(x)$$

で与えられる。sangi の lagrange_interpolation はこの定義式に従い、評価点 $x_i$ 1 つあたり二重ループで $O(n^2)$ の計算量を要する。

Newton 分割差分

同じ多項式を、節点を 1 つずつ追加していく形で書いたのが Newton 形式:

$$p(x) = \sum_{i=0}^{n} c_i \prod_{j=0}^{i-1} (x - x_j)$$

係数 $c_i = f[x_0, \ldots, x_i]$ は分割差分で、次の漸化式で計算する:

$$f[x_i] = y_i, \qquad f[x_i, \ldots, x_{i+k}] = \frac{f[x_{i+1}, \ldots, x_{i+k}] - f[x_i, \ldots, x_{i+k-1}]}{x_{i+k} - x_i}$$

sangi の newton_interpolation_coefficients は、係数配列を $y$ で初期化し、その場で差分を上書きしていく二重ループ ($O(n^2)$) で分割差分表を構成する。 評価は newton_interpolation_evaluate が Horner 法で行い、評価あたり $O(n)$:

$$p(x) = c_0 + (x - x_0)\bigl(c_1 + (x - x_1)\bigl(c_2 + \cdots\bigr)\bigr)$$

Newton 形式は、節点を 1 つ追加するときに既存の係数を再利用できる (差分表を 1 列足すだけ) のが利点。

重心ラグランジュ形式

ラグランジュ補間を、評価あたり $O(n)$ に高速化しつつ数値的に安定化したのが重心 (barycentric) 形式。 まず各節点の重みを

$$w_j = \frac{1}{\displaystyle\prod_{i \neq j} (x_j - x_i)}$$

と定義する (sangi の barycentricLagrangeWeights、$O(n^2)$ で一度だけ計算)。補間多項式は

$$p(x) = \frac{\displaystyle\sum_{j=0}^{n} \frac{w_j}{x - x_j}\, y_j}{\displaystyle\sum_{j=0}^{n} \frac{w_j}{x - x_j}}$$

で表せる (barycentricLagrangeEvaluate)。この第 2 (重心) 形式の利点:

• 評価が $O(n)$: 重みを使い回せるので、同じ節点で多数の点を評価する用途で通常形より速い。
• 数値安定性: 分子と分母に同じ $w_j/(x - x_j)$ が現れ共通因子が約分されるため、桁落ちが起きにくい。 $x$ が節点 $x_j$ に一致すると $1/(x - x_j)$ が発散するが、実装では $x - x_j = 0$ を検出して $y_j$ をそのまま返す。

ただし重心形式も大域多項式である以上、等間隔節点では Runge 現象から逃れられない。 この振動を断ち切るのが、次節以降の区分手法である。

3 次スプライン

各区間 $[x_i, x_{i+1}]$ を 3 次多項式でつなぎ、節点で 1 階・2 階導関数まで連続にしたものが3 次スプライン。 区間 $i$ の多項式を、局所座標 $t = x - x_i$ で

$$S_i(x) = a_i + b_i (x - x_i) + c_i (x - x_i)^2 + d_i (x - x_i)^3$$

と書く。$n$ 区間で係数は $4n$ 個。これを次の条件で決める。

連続条件

• 補間条件: 各区間が両端の節点を通る ($S_i(x_i) = y_i,\ S_i(x_{i+1}) = y_{i+1}$)。
• $C^1$ 連続: 内部節点で 1 階導関数が一致 ($S_{i-1}'(x_i) = S_i'(x_i)$)。
• $C^2$ 連続: 内部節点で 2 階導関数が一致 ($S_{i-1}''(x_i) = S_i''(x_i)$)。

これらだけでは条件が 2 本足りないので、両端で境界条件を課す。sangi の SplineBoundaryCondition は次の 2 種を持つ:

• 自然スプライン (Natural): 両端で $S''(x_0) = S''(x_n) = 0$。端で曲率を 0 にする最も標準的な選択。
• クランプ (Clamped): 両端で 1 階導関数 $S'(x_0), S'(x_n)$ の値を指定する。端の傾きを既知の値に固定したいときに使う。

三重対角系の導出

2 階導関数に対応する係数 $c_i$ ($i = 0, \ldots, n$) を未知数にとると、$C^2$ 連続条件から各内部節点で 1 本ずつ方程式が立つ。 区間幅を $h_i = x_{i+1} - x_i$ として、内部節点 $i$ では

$$h_{i-1}\, c_{i-1} + 2(h_{i-1} + h_i)\, c_i + h_i\, c_{i+1} = 3\!\left(\frac{y_{i+1} - y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}}\right)$$

となる。各式が隣り合う 3 つの未知数 $c_{i-1}, c_i, c_{i+1}$ しか含まないため、係数行列は三重対角になる。境界条件 (自然なら $c_0 = c_n = 0$) を加えると、係数ベクトル $c = (c_0, \ldots, c_n)^\top$ に対する線形系

$$M\,c = d$$

が得られる。$M$ は対角優位な三重対角行列で正則。

Thomas 法による $O(n)$ 求解

三重対角系は、一般の $O(n^3)$ のガウス消去を使わず、Thomas 法 (三重対角専用の前進消去 + 後退代入) で $O(n)$ で解ける。 sangi の cubic_spline_coefficients はこの 2 段で構成される:

• 前進消去: 対角成分を $l_i = 2(h_{i-1} + h_i) - h_{i-1}\,\mu_{i-1}$、上対角を $\mu_i = h_i / l_i$ と更新しつつ、右辺の中間解 $z_i$ を順に求める。
• 後退代入: $c_{n} = z_n$ から逆順に $c_i = z_i - \mu_i\, c_{i+1}$ を解く。あわせて区間係数 $b_i, d_i$ を $b_i = \dfrac{y_{i+1} - y_i}{h_i} - \dfrac{h_i (c_{i+1} + 2 c_i)}{3}$, $d_i = \dfrac{c_{i+1} - c_i}{3 h_i}$ で確定する。

評価は cubic_spline_evaluate が担い、評価点 $x$ が属する区間 $i$ を探してから $S_i(x)$ を Horner 風に計算する。 係数表記 $\{a_i, b_i, c_i, d_i\}$ は他の区分手法 (PCHIP・Akima・Catmull-Rom) と共通なので、これらは同じ評価関数で値を求められる。

形状保存補間

3 次スプラインは最も滑らか ($C^2$) だが、急峻に立ち上がるデータや段差状のデータでは、データにない山や谷 (overshoot) を作ってしまう。 計測値の単調性や符号を崩したくない場面では、滑らかさを少し犠牲にしてでも形状を保存する手法を使う。 以下はいずれも区間ごとの 3 次 Hermite 多項式

$$p_i(t) = y_i\,(1-t)^2(1+2t) + y_{i+1}\,t^2(3-2t) + d_i\,h_i\,t(1-t)^2 + d_{i+1}\,h_i\,t^2(t-1), \quad t = \frac{x - x_i}{h_i}$$

に基づき、各節点の傾き $d_i$ の決め方だけが異なる。傾きを節点間の係数 $\{a,b,c,d\}$ に変換すれば、3 次スプラインと同じ評価関数で扱える。

PCHIP (Fritsch-Carlson)

PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) は、Fritsch-Carlson のアルゴリズムで単調性を保つ手法 (sangi の pchipCoefficients)。 区間傾きを $\delta_i = (y_{i+1} - y_i)/h_i$ として、内部節点の傾きは次のように決める:

• 隣接区間の傾きが同符号 ($\delta_{i-1}\,\delta_i > 0$): 区間幅で重み付けした調和平均 $$d_i = \frac{w_1 + w_2}{\dfrac{w_1}{\delta_{i-1}} + \dfrac{w_2}{\delta_i}}, \quad w_1 = 2h_i + h_{i-1},\ w_2 = h_i + 2h_{i-1}$$ 調和平均なので、片側の傾きが 0 (平坦) なら $d_i = 0$ になる。
• 傾きの符号が異なる、またはいずれかが 0 (極値の節点): $d_i = 0$ とし、その点で山・谷を作らない。

端点は Bessel の片側 3 点公式で初期推定し、$d_0$ と $\delta_0$ の符号が逆なら 0 に、$|d_0|$ が $3|\delta_0|$ を超えるなら $3\delta_0$ に制限する。

最後に全区間で Fritsch-Carlson の単調性条件を課す。 $\alpha = d_i/\delta_i$, $\beta = d_{i+1}/\delta_i$ とおくと、区間が単調になる十分条件は

$$\alpha^2 + \beta^2 \le 9$$

これを破る区間では $\tau = 3/\sqrt{\alpha^2 + \beta^2}$ を掛けて $d_i, d_{i+1}$ を半径 3 の円内に引き戻す。 この縮小により、単調なデータでは補間も単調になり overshoot が起きない。結果は $C^1$ 連続。

修正 Akima

Akima 法 (1970) は、各節点の傾きを近傍 4 区間の局所傾きから、区間傾きの差で重み付けして推定する手法。 外れ値的な 1 点に過敏に反応せず、自然で振動の少ない曲線を与えるのが特徴。 sangi の modifiedAkimaCoefficients は、傾きの重みを

$$w_1 = |\delta_{i+1} - \delta_i| + \varepsilon, \qquad t_i = \frac{w_1\, \delta_{i-1} + w_2\, \delta_i}{w_1 + w_2}$$

のように構成する。元の Akima では分母が 0 になりうる (連続する区間の傾きが等しいとき) 問題を、微小量 $\varepsilon$ を加える修正版で回避し安定化している ($\varepsilon$ は傾き差の最大値に機械精度を掛けたスケール)。 両端は仮想的な区間傾きを線形外挿で補い、内部と同じ式で扱う。結果は $C^1$ 連続。

Catmull-Rom

Catmull-Rom スプラインは、各節点の接線をその両隣の節点から決める $C^1$ 連続な補間で、コンピュータグラフィックスの曲線補間で広く使われる (sangi の catmullRomCoefficients)。 内部節点の傾きは前後の節点を結ぶ割線で与える:

$$t_i = \frac{y_{i+1} - y_{i-1}}{x_{i+1} - x_{i-1}}$$

端点は片側差分を使う。傾き制限を持たないため PCHIP より滑らかに見えるが、単調性は保証しない。 制御点をそのまま通る扱いやすさが利点で、形状を直接デザインする用途に向く。

5 次 Hermite

各節点で値 $f$・1 階導関数 $f'$・2 階導関数 $f''$ までを指定し、区間を 5 次多項式で結ぶのが5 次 (quintic) Hermite 補間 (sangi の quinticHermiteCoefficients)。 区間ごとに 6 係数

$$S_i(t) = a + b\,t + c\,t^2 + d\,t^3 + e\,t^4 + f\,t^5, \quad t = x - x_i$$

を持ち、$C^2$ 連続を達成する。境界条件 $S(0)=y_0,\ S(h)=y_1,\ S'(0)=f'_0,\ \ldots,\ S''(h)=f''_1$ から低次 3 係数 $a, b, c$ が直接定まり、高次 3 係数 $a_3, a_4, a_5$ は $3 \times 3$ の線形系をガウス消去で解いて得る。 導関数情報が手元にあるとき、3 次 Hermite より高い滑らかさと精度が出せる。

なぜ単調データで overshoot しないか

3 次スプラインが overshoot するのは、$C^2$ 連続を満たすために節点の傾き $d_i$ が大域的に連立で決まり、局所的なデータ形状から外れた大きな傾きを取りうるからである。 対して形状保存系は、節点の傾きを局所的に・上限付きで決める:

• 極値の節点 ($\delta_{i-1}$ と $\delta_i$ が逆符号) で傾きを 0 にする → データにない山・谷を作らない。
• 傾きの大きさを区間傾き $\delta_i$ の定数倍以内に制限する (Fritsch-Carlson の $\alpha^2 + \beta^2 \le 9$) → 区間内で 3 次多項式が単調であり続ける。

この「局所性 + 傾き制限」が、単調データに対する単調性 (= overshoot なし) の理論的根拠である。 ただし $C^2$ 連続性は犠牲になり、形状保存系は一般に $C^1$ までしか保証しない。 滑らかさと形状保存はトレードオフの関係にあり、用途に応じて手法を選ぶ。

比較表

大まかな指針: 滑らかさ最優先なら 3 次スプライン、単調性・符号を保ちたいなら PCHIP、外れ値に頑健で自然な曲線が欲しいなら修正 Akima、CG で制御点を直接扱うなら Catmull-Rom、導関数まで既知で高精度が要るなら 5 次 Hermite。 同じ節点で何度も値を引くだけなら重心ラグランジュが速い。

参考文献

• de Boor, C. (2001). A Practical Guide to Splines. Revised ed. Springer.
• Fritsch, F. N., & Carlson, R. E. (1980). "Monotone Piecewise Cubic Interpolation". SIAM Journal on Numerical Analysis, 17(2), 238–246.
• Akima, H. (1970). "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures". Journal of the ACM, 17(4), 589–602.
• Berrut, J.-P., & Trefethen, L. N. (2004). "Barycentric Lagrange Interpolation". SIAM Review, 46(3), 501–517.
• Catmull, E., & Rom, R. (1974). "A Class of Local Interpolating Splines". In Computer Aided Geometric Design (pp. 317–326). Academic Press.