数値微分の実装と精度の限界｜Python

数値微分は、関数の解析的な導関数がわからないとき、あるいは計算が複雑なときに近似値を求める手法だ。アイデアはシンプルだが、刻み幅の選び方によって精度が大きく変わる。浮動小数点数の限界も理解しておく必要がある。

前進差分

微分の定義に基づく最も単純な近似が前進差分だ。

$f^{'} (x) \approx \frac{f ( x + h ) - f ( x )}{h}$

def forward_diff(f, x, h=1e-5):
    """前進差分による数値微分"""
    return (f(x + h) - f(x)) / h

# f(x) = x^2 の微分（解析解は 2x）
f = lambda x: x**2

print(forward_diff(f, 3))  # 6.000009999951316（真の値は 6）

誤差は約 $1 0^{- 5}$ のオーダーで、 $h$ と同程度になる。これは前進差分が一次精度であることを示している。

中心差分

中心差分は前進差分より高精度だ。

$f^{'} (x) \approx \frac{f ( x + h ) - f ( x - h )}{2 h}$

def central_diff(f, x, h=1e-5):
    """中心差分による数値微分"""
    return (f(x + h) - f(x - h)) / (2 * h)

f = lambda x: x**2

print(central_diff(f, 3))  # 5.999999999999339（より正確）

誤差は $h^{2}$ のオーダーになり、二次精度を持つ。同じ $h = 1 0^{- 5}$ でも、前進差分より格段に正確な結果が得られる。

前進差分

一次精度（誤差 $O (h)$ ）。計算が単純だが精度は低い。

中心差分

二次精度（誤差 $O (h^{2})$ ）。計算量はほぼ同じで精度が高い。通常はこちらを使う。

二階微分

二階微分は中心差分を拡張して計算できる。

$f^{''} (x) \approx \frac{f ( x + h ) - 2 f ( x ) + f ( x - h )}{h ^{2}}$

def second_diff(f, x, h=1e-4):
    """二階微分の近似"""
    return (f(x + h) - 2*f(x) + f(x - h)) / (h**2)

# f(x) = x^3 の二階微分（解析解は 6x）
f = lambda x: x**3

print(second_diff(f, 2))  # 12.000000000287557（真の値は 12）

二階微分では $h^{2}$ で割るため、丸め誤差の影響を受けやすい。 $h$ を大きめ（ $1 0^{- 4}$ 程度）に取る必要がある。

刻み幅 h の選び方

$h$ を小さくすれば打ち切り誤差は減るが、丸め誤差が増える。この二つのバランスで最適な $h$ が決まる。

import math

f = lambda x: math.sin(x)
x = 1.0
true_value = math.cos(x)  # sin の微分は cos

for exp in range(-1, -17, -1):
    h = 10 ** exp
    approx = central_diff(f, x, h)
    error = abs(approx - true_value)
    print(f'h = 10^{exp:3d}: error = {error:.2e}')

実行すると、 $h = 1 0^{- 8}$ あたりで誤差が最小になり、それより小さい $h$ では丸め誤差が支配的になって精度が悪化する。

なぜ丸め誤差が問題になるか

$h$ が非常に小さいと、 $f (x + h)$ と $f (x - h)$ の差がごくわずかになる。浮動小数点数では有効桁数が限られているため、この差の計算で精度が失われる（桁落ち）。

import math

f = lambda x: math.sin(x)
x = 1.0

# h が極端に小さい場合
h = 1e-15
print(f(x + h))  # 0.8414709848078965
print(f(x - h))  # 0.8414709848078965
print(f(x + h) - f(x - h))  # 0.0（差がゼロになってしまう）

差がゼロになると微分値もゼロになり、完全に間違った結果になる。

最適な刻み幅の理論値

中心差分の場合、打ち切り誤差と丸め誤差のバランスから、最適な $h$ は次のようになる。

$h_{opt} \approx ϵ^{1/3}$

ここで $ϵ$ は機械イプシロン（約 $2.2 \times 1 0^{- 16}$ ）だ。

import sys

eps = sys.float_info.epsilon
h_opt = eps ** (1/3)
print(f'機械イプシロン: {eps:.2e}')      # 2.22e-16
print(f'最適な h: {h_opt:.2e}')          # 6.06e-06

理論的には $h \approx 1 0^{- 5}$ から $1 0^{- 6}$ あたりが最適ということになる。

高階の差分公式

より高精度な差分公式も存在する。5 点公式は四次精度を持つ。

$f^{'} (x) \approx \frac{- f ( x + 2 h ) + 8 f ( x + h ) - 8 f ( x - h ) + f ( x - 2 h )}{12 h}$

def five_point_diff(f, x, h=1e-3):
    """5 点公式による数値微分（四次精度）"""
    return (-f(x + 2*h) + 8*f(x + h) - 8*f(x - h) + f(x - 2*h)) / (12 * h)

import math
f = lambda x: math.sin(x)
x = 1.0

print(f'中心差分:   {central_diff(f, x, 1e-3):.15f}')
print(f'5点公式:    {five_point_diff(f, x, 1e-3):.15f}')
print(f'真の値:     {math.cos(x):.15f}')

同じ $h$ でも 5 点公式の方が高精度だ。ただし関数の評価回数が増えるため、計算コストとのトレードオフになる。

自動微分との比較

数値微分には本質的に精度の限界がある。高精度が必要な場合や、勾配を大量に計算する機械学習では、自動微分（Automatic Differentiation）が使われる。PyTorch や JAX などのフレームワークが自動微分を提供している。

数値微分

実装が簡単。精度に限界があり、高次元では計算コストが高い。

自動微分

機械精度で正確な勾配が得られる。フレームワークの導入が必要。

数値微分は「とりあえず微分値が欲しい」場面で便利だが、限界を理解した上で使うことが重要だ。