Abstract

みなさん微分していますか？私はしています。数値微分は自動微分のテストに使えるし実装が楽で嬉しいです。なので精度よくやりたくなってきます。ところで前進差分の打ち切り誤差が有限の幅

h

に対して

O (h)

な一方、中心差分は

O (h^{2})

であるというのは割と有名ですが、各

n

に対して

O (h^{n})

で計算する式が常に存在します。(たぶん)
調べたのですが出てこなかったので導出してみます。

前進差分・中心差分の精度の確認

まずは前進差分と中心差分の精度を確認してみます。

前進差分は以下のような計算で得られます。

f

の

x

における

h

差分の前進差分による近似

f_{f o r w a r d}^{'} (x; h)

は、

f_{f o r w a r d}^{'} (x; h) = \frac{f (x + h) - f (x)}{h}

これが一次精度、つまり誤差

e (h)

が

O (h)

であることを確かめます。
上の式の右辺をテイラー展開すると、

\begin{aligned} f_{f o r w a r d}^{'} (x; h) & = \frac{f (x + h) - f (x)}{h} \\ = \frac{f (x) + f^{'} (x) h + \frac{f^{″} (x)}{2!} h^{2} + \frac{f^{‴} (x)}{3!} h^{3} + \dots - f (x)}{h} \\ = \frac{f^{'} (x) h + O (h^{2})}{h} \\ = f^{'} (x) + O (h) \end{aligned}

中心差分についても、全く同様に右辺をテイラー展開することで確かめられます。
Pythonでそれぞれ実装してみます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("module://matplotlib_pyodide.html5_canvas_backend")
def numerical_diff(f, x, h):
    return (f(x+h) - f(x)) / h
def numerical_diff_center(f, x, h):
    return (f(x+h) - f(x-h)) / (2*h)
def f(x):
    return np.sin(x)
h = np.logspace(0, -20, num=100)
grad_forward = numerical_diff(f, np.pi/3, h)
grad_center = numerical_diff_center(f, np.pi/3, h)
grad_true = np.cos(np.pi/3)
error_forward = np.abs(grad_forward - grad_true)
error_center = np.abs(grad_center - grad_true)
plt.plot(h, error_forward, label='forward')
plt.plot(h, error_center, label='center')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('h')
plt.ylabel('grad')
plt.legend()
plt.show()
plt.close()
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

よく言われる、打ち切り誤差と桁落ちのトレードオフが確認できます

任意の k 階 n 次精度数値微分

本題を導出しようと思います。まずは議論のしやすさのために特に

n

が偶数の場合について考えます。

任意の

n, k \leq n + 1

に対して、

n

が偶数のとき、ある

w_{n, k} \in R^{n + 1}

が存在して

p = \frac{n}{2}

として

f_{n}^{(k)} (x; h) = w_{n, k} \cdot {(\begin{matrix} f (x - p h) \\ f (x - (p - 1) h) \\ ⋮ \\ f (x) \\ ⋮ \\ f (x + (p - 1) h) \\ f (x + p h) \end{matrix})}^{⊤}

は

f_{n}^{(k)} (x; h) = f^{(k)} (x) + O (h^{n})

を満たす。

つまり、各

n

k

に対して中心差分と同様に、微分係数を求めたい点から左右に

h

ずつ幅をとって点を取り評価した

n + 1

点に対して、適切な重み付き和が存在して

k

階導関数の

n

次精度の数値微分を実現できるということです。

証明してみます。

f (x + k h) = f (x) + \frac{f^{'} (x)}{1!} k h + \frac{f^{″} (x)}{2!} k^{2} h^{2} + \dots + \frac{f^{(n)} (x)}{n!} k^{n} h^{n} + O (h^{n + 1})

より、

これを

k = - p, - p + 1, \dots, p - 1, p

について足し合わせて

k + 1

番目の項以外の係数を

0

にすればよい。

したがって、

A = (\begin{matrix} 1 & 1 & \dots & 1 & 1 & 1 & \dots & 1 & 1 \\ - p & (- p + 1) & \dots & 1 & 0 & 1 & \dots & (p - 1) & p \\ (- p)^{2} & (- p + 1)^{2} & \dots & (- 1)^{2} & 0 & 1^{2} & \dots & (p - 1)^{2} & p^{2} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ (- p)^{n} & (- p + 1)^{n} & \dots & (- 1)^{n} & 0 & 1^{n} & \dots & (p - 1)^{n} & p^{n} \end{matrix}) \in R^{(n + 1) \times (n + 1)}

と

A

を定めたときに

A w = e_{k + 1}

を満たす

w \in R^{n + 1}

が存在すればよい。
ここで、

A

は ヴァンデルモンド行列 で、

p \neq 0

のとき正則。
したがって、

w = A^{- 1} e_{k + 1}

なる

w

がただ存在して条件を満たす。

正則性について

ヴァンデルモンド行列は、

A = (\begin{matrix} 1 & 1 & \dots & 1 & 1 & 1 & \dots & 1 & 1 \\ x_{1} & x_{2} & \dots & x_{n - 1} & x_{n} & x_{1} & \dots & x_{n - 1} & x_{n} \\ x_{1}^{2} & x_{2}^{2} & \dots & x_{n - 1}^{2} & x_{n}^{2} & x_{1}^{2} & \dots & x_{n - 1}^{2} & x_{n}^{2} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ x_{1}^{n - 1} & x_{2}^{n - 1} & \dots & x_{n - 1}^{n - 1} & x_{n}^{n - 1} & x_{1}^{n - 1} & \dots & x_{n - 1}^{n - 1} & x_{n}^{n - 1} \\ x_{1}^{n} & x_{2}^{n} & \dots & x_{n - 1}^{n} & x_{n}^{n} & x_{1}^{n} & \dots & x_{n - 1}^{n} & x_{n}^{n} \end{matrix}) \in R^{n \times n}

という各列が等比数列のようになっている行列のことを言います。 (各行の派閥もいるらしいです。)
ヴァンデルモンド行列は、その行列式が

det A = \prod_{1 \leq i < j \leq n} (x_{j} - x_{i})

となることが知られていて、

x_{i}

がすべて異なるときには

0

でないことがわかります。
したがって今回の場合は

x_{i} = - p + i

であり、

p \neq 0

のときはすべて異なるので正則であることがわかります。

というわけで確かに存在したうえに一意で、また具体的に計算することができました。
実際に計算してみようと思います。
右から

e_{k}

をかけるのはちょうど

k

列目のベクトルを取り出すことなのでそう実装します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
import numpy as np
import matplotlib.pyplot as plt
def central_weight(n, k):
    p = n // 2
    A = np.zeros((n+1, n+1), dtype=float)
    r = np.arange(-p, p+1)
    
    for i in range(n+1):
        A[i] = r**i
        
    w = np.linalg.inv(A)[:, k]
    return w
def numerical_diff(f, x, n, h=1e-5, k=1):
    w = central_weight(n, k)
    p = n // 2
    
    w = central_weight(n, k)
    x = np.arange(-p, p+1) * h + x
    y = f(x)
    return np.dot(w, y) * np.prod(np.arange(1, k+1)) / h**k
def f(x):
    return np.sin(x)
def df(x):
    return np.cos(x)
x = np.pi / 3
print('w_2,1 =', central_weight(2, 1))
print('w_4,1 =', central_weight(4, 1))
print('w_4,2 =', central_weight(4, 2))
print(numerical_diff(f, x, 2))
print(numerical_diff(f, x, 4))    
print(numerical_diff(f, x, 4, k=2))
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

と、いい感じに計算できています。
桁落ちとのトレードオフも確認してみます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
n = range(2, 12, 2)
x = np.pi / 3
h = np.logspace(0, -10, 100)
label = [f'O(h^{i})' for i in n]
err = np.zeros((len(n), len(h)))
for i, _n in enumerate(n):
    for j, _h in enumerate(h):
        err[i, j] = abs(numerical_diff(f, x, _n, h=_h) - df(x))
plt.figure(figsize=(8, 6))
for i in range(len(n)):
    plt.loglog(h, err[i], label=label[i])
plt.legend()
plt.grid()
plt.xlabel('h')
plt.ylabel('Error')
plt.show()
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

打ち切り誤差の変化によって最適な点が変化していることがわかります。

感想

線形代数が役に立ってすごい
また、上のアルゴリズムは逆行列を求めるところがボトルネックで　時間計算量は

O (N^{3})

　ですが、調べたところによるとヴァンデルモンド行列の逆行列は

O (N^{2})

で求められるようなのでそうすればもう少し早くなりそうです。 (中心差分よりも正確にやろうとする場面すらなかなかみないので実用的ではないですが)

Abstract

前進差分・中心差分の精度の確認

任意の k 階 n 次精度数値微分

感想

今日の一曲