Pythonで学ぶ数値計算のアルゴリズムと実装 第2回：多項式補間とチェビシェフ補間の戦略

Python
2026/01/01
この記事は約22分で読めます。

こんにちは、JS2IIUです。
前回の記事では、既知の複雑な関数を扱いやすい多項式に変換する「テイラー展開」について学びました。しかし、実世界のデータ解析やAI開発においては、元となる関数の形が最初から分かっていることは稀です。多くの場合、私たちが手にしているのは「いくつかの離散的な観測データ（点）」のみです。

これらの点と点の間をどのように埋め、未知の値を推定するか。この課題を解決するのが「補間（Interpolation）」という技術です。今回は、数値計算における補間の基礎から、高次多項式を扱う際に必ず直面する「ルンゲ現象」、そしてその回避策である「チェビシェフ補間」までを、PythonとPyTorchを用いて深く掘り下げていきましょう。今回もよろしくお願いします。

1. はじめに：補間とは何か？

「補間」とは、与えられた複数の点（データ点）をすべて通り、それらの点の間を滑らかにつなぐ関数を求める操作を指します。

ここで注意したいのは、機械学習でよく使われる「回帰（Regression）」との違いです。回帰は、データにノイズが含まれていることを前提とし、必ずしもすべての点を通る必要はありません（全体の傾向を捉えることが目的です）。一方、補間は「与えられたデータ点は正確である」と仮定し、それらをすべて通る曲線を作ります。

画像のリサイズ、欠損データの推定、あるいは物理シミュレーションにおける格子点間の値の取得など、補間技術はあらゆるテックスタックの根底に存在しています。

2. 多項式補間の基礎：ラグランジュ補間

\(n+1\) 個の異なる点 \((x_0, y_0), (x_1, y_1), \dots, (x_n, y_n)\) が与えられたとき、これらすべてを通る \(n\) 次以下の多項式 \(P(x)\) はただ一つ存在します。これを構成する最も直感的な方法が「ラグランジュ補間（Lagrange Interpolation）」です。

ラグランジュ補間多項式 \(P(x)\) は、次の数式で表されます：

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

ここで \(\ell_i(x)\) は「ラグランジュ基本多項式（Lagrange basis polynomial）」と呼ばれ、

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

という形で定義されます。

この \(\ell_i(x)\) の特徴は、\(x = x_i\) のとき \(\ell_i(x_i) = 1\) となり、\(x = x_j\) （\(j \neq i\)）のとき \(\ell_i(x_j) = 0\) となることです。これは、クロネッカーのデルタ \(\delta_{ij}\) の性質と一致します。

$$
\ell_i(x_j) = \delta_{ij} = \begin{cases} 1 & (i = j) \\ 0 & (i \neq j) \end{cases}
$$

この性質により、ラグランジュ補間多項式 \(P(x)\) は与えられたすべての点 \((x_i, y_i)\) を必ず通ります。すなわち、\(P(x_j) = y_j\) となります。

ラグランジュ補間の本質は「各点 \((x_i, y_i)\) に対して、その点でのみ \(1\) となり、それ以外の点では \(0\) となる基底多項式 \(\ell_i(x)\) を構成し、それらに (y_i) を掛けて足し合わせる」ことです。これにより、任意の点での値を滑らかに推定することができます。

ラグランジュ補間のPyTorch実装

以下のコードでは、PyTorchのブロードキャスト機能を活用して、ラグランジュ補間を効率的に計算しています。

Python
import torch
import matplotlib.pyplot as plt

def lagrange_interpolation(x_points, y_points, x_grid):
    """
    ラグランジュ補間を実行する関数

    Args:
        x_points (torch.Tensor): 既知の点のx座標
        y_points (torch.Tensor): 既知の点のy座標
        x_grid (torch.Tensor): 値を推定したいxの範囲（グラフ描画用など）

    Returns:
        torch.Tensor: 推定されたyの値
    """
    n = len(x_points)
    # x_gridの形状に合わせて結果を初期化
    y_pred = torch.zeros_like(x_grid)

    for i in range(n):
        # ラグランジュ基底多項式 L_i(x) の計算
        li = torch.ones_like(x_grid)
        for j in range(n):
            if i != j:
                # 自身の点以外では0になるように積を計算
                li *= (x_grid - x_points[j]) / (x_points[i] - x_points[j])

        # y_i * L_i(x) を加算
        y_pred += y_points[i] * li

    return y_pred

# テスト用のデータ点
x_obs = torch.tensor([-1.0, 0.0, 1.0, 2.0])
y_obs = torch.tensor([1.0, 0.0, 1.0, 8.0]) # y = x^3 の一部のようなデータ

# 描画用の細かいグリッド
x_range = torch.linspace(-1.2, 2.2, 100)
y_interp = lagrange_interpolation(x_obs, y_obs, x_range)

# 可視化
plt.figure(figsize=(8, 5))
plt.plot(x_range, y_interp, label='Lagrange Interpolation', color='blue')
plt.scatter(x_obs, y_obs, color='red', label='Observed Points')
plt.legend()
plt.title('Basic Lagrange Interpolation')
plt.grid(True)
plt.show()

この実装は理解しやすいですが、計算量は \(O(n^2)\) であり、新しい点が追加されるたびにすべての基底多項式を計算し直さなければならないという欠点があります。

3. より効率的なアプローチ：ニュートン補間

計算効率の面でラグランジュ補間を上回るのが「ニュートン補間（Newton Interpolation）」です。この手法では「差分商（Divided Differences）」と呼ばれる値を計算し、多項式を積み上げていく形式をとります。

ニュートン補間の最大の特徴は、新しいデータ点が一つ追加されたとき、既存の計算結果を無駄にせず、新しい項を追加するだけで済む点にあります。これはリアルタイムでデータが増えていくストリーミングデータの処理などにおいて非常に有利です。

ニュートン補間の実装：差分商表の作成

ニュートン補間の係数を求めるために、まずは差分商表を作成する必要があります。

Python
def compute_divided_diffs(x_points, y_points):
    """
    ニュートン補間に必要な差分商表を計算する
    """
    n = len(x_points)
    # 差分商テーブルの初期化
    coef = y_points.clone().detach().float()

    for j in range(1, n):
        # 下から順に更新することで、1次元配列でメモリを節約しながら計算可能
        for i in range(n - 1, j - 1, -1):
            coef[i] = (coef[i] - coef[i-1]) / (x_points[i] - x_points[i-j])

    return coef

def newton_interpolation(coef, x_points, x_grid):
    """
    差分商係数を用いて、特定のxにおける値を計算する
    ホーナー法に似た形式で効率的に計算
    """
    n = len(coef)
    y_pred = torch.ones_like(x_grid) * coef[-1]

    # ホーナー法的な入れ子構造で多項式を評価
    for i in range(n - 2, -1, -1):
        y_pred = y_pred * (x_grid - x_points[i]) + coef[i]

    return y_pred

# 実行例
divided_diffs = compute_divided_diffs(x_obs, y_obs)
y_newton = newton_interpolation(divided_diffs, x_obs, x_range)

ここで、前回の「ホーナー法」の考え方が再登場していることに注目してください。数値計算では、このように過去に学んだ「効率化のパターン」が別のアルゴリズムの中で形を変えて何度も現れます。

Python
import matplotlib.pyplot as plt

# 可視化
plt.figure(figsize=(8, 5))
plt.plot(x_range, y_newton, label='Newton Interpolation', color='green')
plt.scatter(x_obs, y_obs, color='red', label='Observed Points')
plt.legend()
plt.title('Newton Interpolation')
plt.grid(True)
plt.savefig('newton.png')
plt.show()

4. 高次補間の罠：ルンゲ現象の視覚化

さて、ここまで「多項式の次数を上げれば、より正確にデータを再現できる」と考えてきたかもしれません。しかし、ここには数値計算における非常に有名な罠が潜んでいます。それが「ルンゲ現象（Runge’s Phenomenon）」です。

ルンゲ現象とは、高次の多項式で等間隔な点を補間しようとすると、区間の端点付近で関数の値が激しく振動してしまう現象を指します。

これを体験するために、有名なルンゲ関数 \(f(x) = \frac{1}{1 + 25x^2}\) を使って実験してみましょう。

Python
def runge_function(x):
    return 1 / (1 + 25 * x**2)

# 次数を11（データ点12個）に設定
n_points = 12
x_equi = torch.linspace(-1, 1, n_points)
y_equi = runge_function(x_equi)

x_fine = torch.linspace(-1, 1, 200)
y_true = runge_function(x_fine)
y_interp_runge = lagrange_interpolation(x_equi, y_equi, x_fine)

plt.figure(figsize=(10, 6))
plt.plot(x_fine, y_true, 'k--', label='True Runge Function')
plt.plot(x_fine, y_interp_runge, 'r', label='Equidistant Interpolation (n=11)')
plt.scatter(x_equi, y_equi, color='red')
plt.ylim(-0.5, 1.5)
plt.title("Runge's Phenomenon: Oscillations at the edges")
plt.legend()
plt.show()

グラフを見ると、中央付近はうまく近似できていますが、両端（\(x = -1, 1\) 付近）では真の関数から大きく外れ、激しく上下に振れていることがわかります。これでは実用的な予測には使えません。なぜこのようなことが起きるのでしょうか。

直感的に説明すると、等間隔に点を配置した場合、多項式の重みが端の点に対して過剰に反応してしまうためです。これを解決するには、「点の配置（サンプリング戦略）」を見直す必要があります。

5. 賢いサンプリング：チェビシェフ補間

ルンゲ現象を回避し、補間誤差を最小化するための魔法の杖が「チェビシェフノード（Chebyshev Nodes）」です。

チェビシェフノードの戦略は、「誤差が出やすい端点付近に、より多くのサンプリング点を配置する」というものです。数学的には、半円を等間隔な角度で区切り、その点をx軸に投影した位置に対応します。これにより、多項式が端で暴走するのを抑え込むことができます。

\([-1, 1]\) の区間における \(n\) 個のチェビシェフノードは、以下の式で求められます。

$$x_i = \cos\left( \frac{2i + 1}{2n} \pi \right), \quad i = 0, \dots, n-1$$

チェビシェフ補間の実装と比較

それでは、等間隔サンプリングとチェビシェフノードサンプリングを比較してみましょう。

Python
def get_chebyshev_nodes(n, a=-1, b=1):
    """
    区間[a, b]におけるn個のチェビシェフノードを生成する
    """
    i = torch.arange(n).float()
    # 基本のチェビシェフノード（[-1, 1]）
    nodes = torch.cos((2 * i + 1) / (2 * n) * torch.pi)
    # 任意の区間[a, b]にスケーリング
    return 0.5 * (a + b) + 0.5 * (b - a) * nodes

# チェビシェフノードでのサンプリング
x_cheb = get_chebyshev_nodes(n_points)
y_cheb = runge_function(x_cheb)

# 補間の実行
y_interp_cheb = lagrange_interpolation(x_cheb, y_cheb, x_fine)

# 比較の可視化
plt.figure(figsize=(10, 6))
plt.plot(x_fine, y_true, 'k--', label='True Runge Function')
plt.plot(x_fine, y_interp_runge, 'r:', label='Equidistant (Runge Phenomenon)')
plt.plot(x_fine, y_interp_cheb, 'g', label='Chebyshev Nodes (Stable)', lw=2)
plt.scatter(x_cheb, y_cheb, color='green', marker='x')
plt.ylim(-0.2, 1.2)
plt.title("Comparison: Equidistant vs Chebyshev Nodes")
plt.legend()
plt.show()

結果は一目瞭然です。チェビシェフノードを用いると、同じ11次の多項式であっても、端点での振動が劇的に抑えられ、真の関数を驚くほど正確にトレースできています。

6. 実践的なアドバイス：SciPyの活用

教育的な目的でアルゴリズムを自作してきましたが、実際のプロダクション環境では、より洗練されたライブラリを使用することが推奨されます。

Pythonの科学計算ライブラリ SciPy には、scipy.interpolate.BarycentricInterpolator というクラスがあります。これは今回学んだラグランジュ補間の改良版（重心ラグランジュ補間）を実装しており、数値的な安定性が非常に高く、計算量も抑えられています。

Python
from scipy.interpolate import BarycentricInterpolator

# SciPy版の実装例（NumPy配列を使用）
interp_fn = BarycentricInterpolator(x_cheb.numpy(), y_cheb.numpy())
y_scipy = interp_fn(x_fine.numpy())

自作の実装で原理を理解し、実務ではこうした最適化されたライブラリを使い分けるのが、エンジニアとしての賢い選択です。

まとめ

第2回となる今回は、離散データの間を埋める「補間」について学びました。

  • ラグランジュ補間は概念として分かりやすく、基底多項式の足し合わせで構成される。
  • ニュートン補間は差分商表を用いることで、データの追加に対して効率的に対応できる。
  • ルンゲ現象は、等間隔サンプリングによる高次多項式補間が引き起こす致命的な誤差である。
  • チェビシェフ補間は、ノードの配置を工夫することでルンゲ現象を回避し、最良に近い近似を実現する。

データの間を「単につなぐ」だけではなく、「どうつなげば誤差が少ないか」という視点を持つことは、より精度の高いAIモデルやシミュレータを構築する上で欠かせません。

次回は、これらの関数近似技術をさらに応用し、面積を求める「数値積分」の世界へ足を踏み入れます。中点則、台形則といった基本から、現代的な高度な積分公式までを網羅する予定です。お楽しみに。

最後まで読んでいただきありがとうございました。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装 第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装 第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装 第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装 第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装 第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装 第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装 第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装 第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装 第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装 第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装 第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装 第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/

コメント

タイトルとURLをコピーしました