あざらしとペンギンの問題

主に漫画、数値計算、幾何計算、TCS、一鰭旅、水族館、鰭脚類のことを書きます。

アクセルほど微分精度は加速していくよ!

今回のテーマは微分の計算です。やりたいことは、実数の範囲で定義されたなめらかな関数 f(x) と実数値 a が与えられたとき、fx = a での微分係数 f'(a) を求めることです。

もし f が初等関数の式として与えられているならば話は簡単です。f微分公式を適用していけば必ず導関数 f' が初等関数の式として得られるので、それに x = a を代入して計算すればいいというだけです。もちろんこの計算はコンピュータで実行できます。例えば、MathematicaフリーソフトMaxima のような数式処理システムを使えば導関数の計算から値の計算までをやってくれます。さらには、ソフトウェアをインストールしなくても、Web サービスの Wolfram Alpha を利用することができます。実際 Wolfram alpha は大変良いサービスで、大雑把な入力に対しても大抵欲しい答えを返してくれます。

では f の式が与えられていない場合はどうやって計算すればいいか?精度を高めるにはどうすれば?というのが今回の話の主題です。

今回の話では、関数 f(x) := \sin 2xx = x_0 := \pi/8 での微分係数を求めるとして試してみましょう。導関数は公式によって f'(x) = 2 \cos 2x と求められるので、求めたい値の答えは d_0 := f'(\pi/8) = \sqrt{2} = 1.4142135623730950488...\dots です。

単純な方法

定義式をそのまま使う

まず微分の定義からおさらいしましょう。x_0 の近傍で定義された関数 f について、極限

\displaystyle \lim_{h \to 0} \frac{f(x0 + h) - f(a)}{h}

が存在するならば、その極限値x_0 での微分係数といい、f'(x_0) で表すのでした。

というわけで f'(x_0) を計算する最も単純な方法は、小さい h を代入して

\displaystyle D_\mathrm{f}(x, h) := \frac{f(x + h) - f(x)}{h}

を計算することでしょう。

hh = 1 / 2^3 = 0.125 から初めて 1/2 ずつにしていくとし、それぞれの h について D_\mathrm{f}(x_0, h) を計算します。

この記事は結構何度も書き直していて、いくつかのツールを取っ替え引っ替えしてきましたが、今度はごまかし無しにやるためにC言語を使うことにします。数値計算向けのツールについては、過去に使っていたものを含め、記事の最後にいくつか挙げておきました。コードの実行はオンラインの処理系である ideone.com で行うことにします。

Ideone.com - Online Compiler and IDE >> C/C++, Java, PHP, Python, Perl and 40+ other compilers and interpreters

#include <stdio.h>
#include <math.h>

// 計算する値の個数
#define N 40

// 円周率
#define M_PI acos(-1.0)

// 微分係数を計算する点
#define x0 (M_PI / 8)

// 真の値 = sqrt(2)
#define d0 sqrt(2)

// 関数 f
double f(double x) { return sin(2 * x); }

// 相対誤差を計算する関数
double err(double x) { return fabs(x - d0) / d0; }

// D_f
double dif_f(double x, double h) { return (f(x + h) - f(x)) / h; }

int main(void) {
    for(int k = 0; k < N; k++) {
        // 刻み幅
        double h = pow(2, -3 - k);

        // D_f の計算値
        double y = dif_f(x0, h);

        // 値と誤差を表示
        printf("%2d: %e  %.16lf  %e\n", k, h, y, err(y));
    }

    return 0;
}

これを ideone.com で実行した結果は次のようになりました。各列は左から順に hD_\mathrm{f}(h)、真の値との相対誤差 |D_\mathrm{f}(x_0, h) - d_0| / d_0 です。

コードと実行結果

 0: 1.250000e-01  1.2236702388976219  1.347345e-01
 1: 6.250000e-02  1.3222602721296628  6.502080e-02
 2: 3.125000e-02  1.3691132406553272  3.189074e-02
 3: 1.562500e-02  1.3918881069182518  1.578648e-02
 4: 7.812500e-03  1.4031076999036998  7.853031e-03
 5: 3.906250e-03  1.4086749326625352  3.916403e-03
 6: 1.953125e-03  1.4114478334930141  1.955666e-03
 7: 9.765625e-04  1.4128315957475479  9.771980e-04
 8: 4.882812e-04  1.4135228036789158  4.884402e-04
 9: 2.441406e-04  1.4138682392012925  2.441804e-04
10: 1.220703e-04  1.4140409148340041  1.220802e-04
11: 6.103516e-05  1.4141272421165922  6.103764e-05
12: 3.051758e-05  1.4141704031244444  3.051820e-05
13: 1.525879e-05  1.4141919829708058  1.525894e-05
14: 7.629395e-06  1.4142027727357345  7.629426e-06
15: 3.814697e-06  1.4142081675818190  3.814693e-06
16: 1.907349e-06  1.4142108649830334  1.907343e-06
17: 9.536743e-07  1.4142122137127444  9.536469e-07
18: 4.768372e-07  1.4142128881067038  4.767783e-07
19: 2.384186e-07  1.4142132252454758  2.383852e-07
20: 1.192093e-07  1.4142133938148618  1.191887e-07
21: 5.960464e-08  1.4142134785652161  5.926112e-08
22: 2.980232e-08  1.4142135232686996  2.765098e-08
23: 1.490116e-08  1.4142135456204414  1.184592e-08
24: 7.450581e-09  1.4142135530710220  6.577559e-09
25: 3.725290e-09  1.4142135679721832  3.959153e-09
26: 1.862645e-09  1.4142135977745056  2.503258e-08
27: 9.313226e-10  1.4142136573791504  6.717943e-08
28: 4.656613e-10  1.4142136573791504  6.717943e-08
29: 2.328306e-10  1.4142136573791504  6.717943e-08
30: 1.164153e-10  1.4142141342163086  4.043542e-07
31: 5.820766e-11  1.4142150878906250  1.078704e-06
32: 2.910383e-11  1.4142150878906250  1.078704e-06
33: 1.455192e-11  1.4142150878906250  1.078704e-06
34: 7.275958e-12  1.4142150878906250  1.078704e-06
35: 3.637979e-12  1.4142150878906250  1.078704e-06
36: 1.818989e-12  1.4142456054687500  2.265789e-05
37: 9.094947e-13  1.4143066406250000  6.581626e-05
38: 4.547474e-13  1.4143066406250000  6.581626e-05
39: 2.273737e-13  1.4145507812500000  2.384498e-04

上から見ると h が小さくなるほど誤差が小さくなっているのが分かります。しかし、下の方はむしろ誤差が大きくなっているのが分かるでしょうか?特に、下の方は末尾に0が続いていて、明らかに何らかのトラブルが起こっていることを示唆します。

これは『桁落ち』という現象によるもので、精度を高めようと h をあまり小さくし過ぎると、近接する数の引き算によって有効数字のキャンセルが起こり、かえって誤差が大きくなってしまうというジレンマがあるのです。そのため、上の表では小数第7位までが意味を持つ程度の値しか得られていないのです。double が IEEE 754 の float64 として実装されているならば*1、10進の有効桁数はおよそ16桁です。それと比較すると半分くらいの桁数の精度の値しか得られていないということになります。

桁落ちと関係する現象に『情報落ち』というものがあります。これは、計算機上での有限桁の数の計算において情報が失われるというものです。それが明確に現れているのが下の方の値です。最後の値では12桁分しか情報が保たれていないことが判ります。ここまで来るとどう頑張ってもそれ以上精度を高めることはできません。

とにもかくにも一度落ちてしまった精度を取り戻すことはできません。

というわけで、微分を精度よく計算するには、h がそれほど小さくないうち、つまり、桁落ちが酷くならないうちから、真の値に素早く近づくようにすることが必要であることが分かりました。

中心差分近似

微分の精度を高める方法として古くから知られているものに中心差分近似があります。これは、上の D_\mathrm{f}(x, h) の代わりに

\displaystyle D_\mathrm{m}(x, h) := \frac{f(x + h) - f(x  - h)}{2h}

を計算するというものです。こうする方が良さそうだということは、直感的には次の図から見て取れるでしょう。

f:id:azapen6:20151222103901p:plain

先ほどのプログラムの関数 dif_f を次のものに置き換えて計算します。

// D_m
double dif_m(double x, double h) { return (f(x + h) - f(x - h)) / (2 * h); }

この式を使って N = 20 までの範囲で計算した結果は次のようになりました。各列は左から順に hD_\mathrm{m}(x_0, h)、真の値との相対誤差 |D_\mathrm{m}(x_0, h) - \sqrt{2}| / \sqrt{2} で、比較のために4列目に先程の |D_\mathrm{f}(x0, h) - \sqrt{2}| / \sqrt{2} を記しました。

コードと実行結果

 0: 1.250000e-01  1.3995281382501878  1.038416e-02  1.347345e-01
 1: 6.250000e-02  1.4105335907091092  2.602133e-03  6.502080e-02
 2: 3.125000e-02  1.4132930302282549  6.509145e-04  3.189074e-02
 3: 1.562500e-02  1.4139833956233261  1.627525e-04  1.578648e-02
 4: 7.812500e-03  1.4141560185783746  4.068961e-05  7.853031e-03
 5: 3.906250e-03  1.4141991762926978  1.017250e-05  3.916403e-03
 6: 1.953125e-03  1.4142099658447762  2.543130e-06  1.955666e-03
 7: 9.765625e-04  1.4142126632405052  6.357828e-07  9.771980e-04
 8: 4.882812e-04  1.4142133375898993  1.589457e-07  4.884402e-04
 9: 2.441406e-04  1.4142135061772478  3.973646e-08  2.441804e-04
10: 1.220703e-04  1.4142135483243692  9.933949e-09  1.220802e-04
11: 6.103516e-05  1.4142135588608653  2.483522e-09  6.103764e-05
12: 3.051758e-05  1.4142135614947620  6.210754e-10  3.051820e-05
13: 1.525879e-05  1.4142135621550551  1.541776e-10  1.525894e-05
14: 7.629395e-06  1.4142135623187642  3.841781e-11  7.629426e-06
15: 3.814697e-06  1.4142135623696959  2.403655e-12  3.814693e-06
16: 1.907349e-06  1.4142135623551439  1.269341e-11  1.907343e-06
17: 9.536743e-07  1.4142135623842478  7.886103e-12  9.536469e-07
18: 4.768372e-07  1.4142135623842478  7.886103e-12  4.767783e-07
19: 2.384186e-07  1.4142135623842478  7.886103e-12  2.383852e-07

先ほどより早く真の値に近づいていることが見て取れます。精度も8桁程度から12桁程度まで改善されました。上の方を見ると誤差の指数の減り方が D_\mathrm{m}(x_0, h) の方が2倍ほど速くなっています。下の方では桁落ちによって誤差が大きくなっていますが、それでも D_\mathrm{m}(x_0, h) の方が誤差が小さくなっています。

では、どうして中心差分近似の方が定義式にそのまま入れるより精度がよくなるのでしょうか?それを見るには少しばかり計算が必要です。

仮定として f は3回微分可能だとします。このとき、テイラーの定理より定数 A, B が存在して、h \to 0

\displaystyle f(x_0 + h) = f(x_0) + f'(x_0) h + A h^2 + B h^3 + o(h^3)

が成り立ちます。これより

\displaystyle D_\mathrm{f}(x_0, h) = f'(x_0) + A h + B h^2 + o(h^2)

となります。すなわち、誤差のオーダーは h の1次です。

一方、中心差分近似では、

\displaystyle f(x_0 + h) = f(x_0) + f'(x_0) h + A h^2 + B h^3 + o(h^3)
\displaystyle f(x_0 - h) = f(x_0) - f'(x_0) h + A h^2 - B h^3 + o(h^3)

の差を取ると

\displaystyle f(x_0 + h) - f(x_0 - h) = 2 f'(x_0) h + 2 B h^3 + o(h^3)

となることより

\displaystyle D_\mathrm{m}(x_0, h) = f'(x_0) + B h^2 + o(h^2)

となり、誤差のオーダーは h の2次です。よって、こちらの方が h を小さくしていくと誤差が早く0に近づきます。

これで中心差分近似が定義式をそのまま計算するより理論的にも優れた方法であることが分かりました。

加速法

今までは前フリで本題はここからです。以上の単純な方法よりも高い精度を得たいときどうするかを考えましょう。

一般的な設定から導出

まず、一般的な設定として、関数 A(h)h \to 0

\displaystyle A(h) = A + a_1 h^{p_1} + a_2 h^{p_2} + \dots + a_m h^{p_m} + o(h^{p_m})
ただし p_1 < p_2 < \dots < p_m

という展開を持つとします。このとき、極限 \lim_{h \to 0} A(h) = A を求めたいとします。(ただし、h = 0 で評価することはできないとします。)

先ほどと同じように h を半分ずつにしていくことを考えましょう。このとき、

\displaystyle A(h / 2) = A + \frac{a_1 h^{p_1}}{2^{p_1}} +  \frac{a_2 h^{p_2}}{2^{p_2}} + \dots + \frac{a_m h^{p_m}}{2^{p_m}} + o(h^{p_m})

となります。ここで次の計算をすると h^{p_1} の項を消すことができます。

\displaystyle 2^{p_1} A(h / 2) -  A(h) = (2^{p_1} - 1) A + a'_2 h^{p_2} + \dots + a'_m h^{p_m} + o(h^{p_m})
ただし a'_i \; (i=2, 3, \dots, m) はこの計算における h^{p_i} の係数

さらに A を求めるために変形すると、

\displaystyle A^{(1)}(h) := \frac{2^{p_1} A(h / 2) - A(h)}{2^{p_1} - 1}
\displaystyle = A + a''_2 h^{p_2} + \dots + a''_m h^{p_m} + o(h^{p_m})
ただし a''_i \; (i=2, 3, \dots, m) はこの計算における h^{p_i} の係数

となります。上で定義される関数 A^{(1)}(h)A(h) と同じく h \to 0A に収束します。ここで A(h) では誤差が hp_1 次だったのに対し、A^{(1)}(h) では hp_2 次となるため、より速く A に収束することが期待されます。

このようにして収束を早める方法を加速法といいます。ここで説明した方法はリチャードソン補外法(Richardson extrapolation)と呼ばれています。名前は20世紀初期にこの方法を発明したイギリスの数学者リチャードソン*2から取られていますが、それより200年ほど前に和算家の建部賢弘*3が同様の方法を発明していたことが知られています。

加速法は繰り返し適用することによってさらに収束を早めることができます。A^{(1)}(h) の係数を改めて

\displaystyle A^{(1)}(h) = A + a_2 h^{p_2} + a_3 h^{p_3} + \dots + a_m h^{p_m} + o(h^{p_m})

とします。ここで先と同じように

\displaystyle A^{(2)}(h) := \frac{2^{p_2} A^{(1)}(h / 2) - A^{(1)}(h)}{2^{p_2} - 1}

を定義します。すると

\displaystyle A^{(2)}(h) = A + a'_3 h^{p_3} + \dots + a'_m h^{p_m} + o(h^{p_m})
ただし a'_i \; (i=3, 4, \dots, m) はこの計算における h^{p_i} の係数

となります。誤差は hp_3 次になりました。

さらに

\displaystyle A^{(3)}(h) := \frac{2^{p_3} A^{(2)}(h / 2) - A^{(2)}(h)}{2^{p_3} - 1}

を計算すると、誤差のオーダーはさらに小さくなると考えられます。

計算方法

それでは A^{(j)}(h) を計算する方法を考えましょう。記号の統一のため A^{(0)}(h) := A(h) とします。h は適当に小さすぎない h_0 (先の例では h_0 = 2^{-3} = 0.125)から出発して半分ずつにしていくとします。すなわち h_{k + 1} = h_k / 2 です。このとき、

\displaystyle a^{(j)}_k := A^{(j)}(h_k)

と定義すると、

\displaystyle a^{(j + 1)}_k = A^{(j + 1)}(h_k) = \frac{2^{p_j} A^{(j)}(h_k / 2) - A^{(j)}(h)}{2^{p_j} - 1}
\displaystyle = \frac{2^{p_j} A^{(j)}(h_{k + 1}) - A^{(j)}(h)}{2^{p_j} - 1} = \frac{2^{p_j} a^{(j)}_{k+1} - a^{(j)}_k}{2^{p_j} - 1}

すなわち

\displaystyle a^{(j + 1)}_k = \frac{2^{p_j} a^{(j)}_{k+1} - a^{(j)}_k}{2^{p_j} - 1}

となり、a^{(j + 1)}_k を計算するには a^{(j)}_{k+1}, a^{(j)}_k, p_j があればいいということになります。

a^{(j)}_k \; (k = 1, 2, \dots, n) のリストと p = p_j を受け取って a^{(j + 1)}_k \; (k = 1, 2, \dots, n - 1) のリストを計算する関数を次に示します。

// 結果は引数として受け取った配列 b に格納する
// 配列 a の長さは L で指定され、b[0] から b[L-2] までに計算値を入れる
void accel(double a[], int L, int p, double b[]) {
    for (int k = 0; k < L - 1; k++)
        b[i] = (pow(2, p) * a[k + 1] - a[k]) / (pow(2, p) - 1)
}

未知となっている p_j を求めるには、次の計算をします。

\displaystyle \frac{a^{(j)}_{k + 2} - a^{(j)}_{k + 1}}{a^{(j)}_{k + 1} - a^{(j)}_k} = \frac{(2^{-2p_j} - 2^{-p_j}) a_j h^{p_j} + o(h^{p_j})}{(2^{-p_j} - 1) a_j h^{p^j} + o(h^{p_j})}
\displaystyle \approx \frac{2^{-2p_j} - 2^{-p_j}}{2^{-p_j} - 1} = 2^{-p_j}

これより

\displaystyle -\log_2 \frac{a^{(j)}_{k + 2} - a^{(j)}_{k + 1}}{a^{(j)}_{k + 1} - a^{(j)}_k}

を最も近い整数に丸めたものを p_j として採用すればよいでしょう。上の値を表示する関数 order を次に示します。計算値を求めて返すのではなく、敢えて表示するようにしたのは、この計算においてはクリティカルな算術エラー(オーバーフロー、0除算)が発生するためです。

void order(double a[], double L) {
    for (int k = 0; k < L - 2; k++)
        printf("%2d: %lf\n", k, -log2((a[k + 2] - a[k + 1]) / (a[k + 1] - a[k])))
}

微分を加速する

それでは本題の微分を加速しましょう。目標は f(x) := sin 2x について f'(1) を求めることです。加速する元の数列は a^{(0)}_k = D_\mathrm{m}(h_k)h_k = 2^{-k - 3} です。

中心差分近似の議論と同様にして、f2m + 1微分可能なとき

\displaystyle D_\mathrm{m}(h) = f'(a) + a_2 h^2 + a_4 h^4 + \dots + a_{2m} h^{2m} + o(h^{2m})

となります。これより、p_1 = 2 となることが予想されます。これを先の関数 order で確かめてみましょう。

コードと実行結果

 0: 1.995771
 1: 1.998943
 2: 1.999736
 3: 1.999934
 4: 1.999983
 5: 1.999996
 6: 1.999999
 7: 2.000000
 8: 1.999990
 9: 2.000039
10: 2.000125
11: 1.996020
12: 2.011973
13: 1.684498
14: -nan
15: -nan
16: inf
17: nan

これより、最初の方ではちゃんと期待通りの値が出ていることが分かります。よって、まずは p_1 = 2 として accel を呼び、a^{(1)}_k のリストを計算します。

コードと実行結果

 0: 1.250000e-01  1.4142020748620829  8.122897e-06
 1: 6.250000e-02  1.4142128434013035  5.083898e-07
 2: 3.125000e-02  1.4142135174216832  3.178545e-08
 3: 1.562500e-02  1.4142135595633907  1.986761e-09
 4: 7.812500e-03  1.4142135621974721  1.241843e-10
 5: 3.906250e-03  1.4142135623621357  7.749505e-12
 6: 1.953125e-03  1.4142135623724148  4.810763e-13
 7: 9.765625e-04  1.4142135623730308  4.553268e-14
 8: 4.882812e-04  1.4142135623730308  4.553268e-14
 9: 2.441406e-04  1.4142135623734096  2.223251e-13
10: 1.220703e-04  1.4142135623730308  4.553268e-14
11: 6.103516e-05  1.4142135623727274  2.600073e-13
12: 3.051758e-05  1.4142135623751528  1.455005e-12
13: 1.525879e-05  1.4142135623733338  1.687849e-13
14: 7.629395e-06  1.4142135623866732  9.601115e-12
15: 3.814697e-06  1.4142135623502934  1.612328e-11
16: 1.907349e-06  1.4142135623939491  1.474599e-11
17: 9.536743e-07  1.4142135623842478  7.886103e-12
18: 4.768372e-07  1.4142135623842478  7.886103e-12

3列目の誤差を見ると、この時点でかなり速く真の値に近づいていることが分かります。

今選んだ p_1 が適切であることを見るために、p_1 = 1 に変えて order を試してみます。

コードと実行結果

 0: 1.250000e-01  1.4215390431680306  5.179897e-03
 1: 6.250000e-02  1.4160524697474006  1.300304e-03
 2: 3.125000e-02  1.4146737610183973  3.254096e-04
 3: 1.562500e-02  1.4143286415334231  8.137325e-05
 4: 7.812500e-03  1.4142423340070209  2.034462e-05
 5: 3.906250e-03  1.4142207553968547  5.086236e-06
 6: 1.953125e-03  1.4142153606362342  1.271564e-06
 7: 9.765625e-04  1.4142140119392934  3.178913e-07
 8: 4.882812e-04  1.4142136747645964  7.947279e-08
 9: 2.441406e-04  1.4142135904714905  1.986857e-08
10: 1.220703e-04  1.4142135693973614  4.966906e-09
11: 6.103516e-05  1.4142135641286586  1.241371e-09
12: 3.051758e-05  1.4142135628153483  3.127202e-10
13: 1.525879e-05  1.4142135624824732  7.734197e-11
14: 7.629395e-06  1.4142135624206276  3.361050e-11
15: 3.814697e-06  1.4142135623405920  2.298317e-11
16: 1.907349e-06  1.4142135624133516  2.846562e-11
17: 9.536743e-07  1.4142135623842478  7.886103e-12
18: 4.768372e-07  1.4142135623842478  7.886103e-12

誤差を見ると今ひとつ収束が悪いのが分かるでしょう。

話を戻して、先程求めた a[1] に対して再び order を計算し、p_2 の値を決めます。理論的には p_2 = 4 となるはずです。

コードと実行結果

 0: 3.997887
 1: 3.999472
 2: 3.999877
 3: 3.999706
 4: 4.001734
 5: 4.060754
 6: inf

早々とエラーになったのは誤差が非常に近接しているためです。ともかく p_2 = 4 が決まりました。先ほどと同じように accel を呼び、a^{(2)}_k のリストを計算します。以下では見やすさのため最初の10個だけ表示します。

コードと実行結果

 0: 1.250000e-01  1.4142135613039182  7.560223e-10
 1: 6.250000e-02  1.4142135623563752  1.182280e-11
 2: 3.125000e-02  1.4142135623728378  1.819737e-13
 3: 1.562500e-02  1.4142135623730776  1.240373e-14
 4: 7.812500e-03  1.4142135623731134  1.287476e-14
 5: 3.906250e-03  1.4142135623731000  3.454203e-15
 6: 1.953125e-03  1.4142135623730718  1.648597e-14
 7: 9.765625e-04  1.4142135623730308  4.553268e-14
 8: 4.882812e-04  1.4142135623734347  2.400671e-13
 9: 2.441406e-04  1.4142135623730057  6.327473e-14

2番目の時点で既に中心差分近似で得られた最良の値よりも誤差が小さくなっています。さらに、4番目は64ビット浮動小数点型でで到達できる精度の限界にかなり迫っています。一方で、加速するほど桁落ちのために列を先に進めるとすぐ精度が落ちてしまうという問題も見えてきます。

さらに加速を試みて order を計算すると、

コードと実行結果

 0: 5.998424
 1: 6.101168
 2: 2.745899
 3: -nan
 4: -1.081794
 5: -0.542697
 6: -nan

となり、理論上は p_3 = 6 となるところが結果が安定しません。というのも、a[2] の収束が非常に速いため、早い段階で著しい桁落ちが生じてしまうのです。とはいえ最初の二つの値の差は 0.1 程度しかないので、それらを丸めて p_3 = 6 としてもう一度くらい加速してみましょう。

コードと実行結果

 0: 1.250000e-01  1.4142135623730807  1.020560e-14
 1: 6.250000e-02  1.4142135623730989  2.669157e-15
 2: 3.125000e-02  1.4142135623730814  9.734573e-15
 3: 1.562500e-02  1.4142135623731140  1.334579e-14
 4: 7.812500e-03  1.4142135623731000  3.454203e-15
 5: 3.906250e-03  1.4142135623730714  1.679999e-14
 6: 1.953125e-03  1.4142135623730301  4.600371e-14
 7: 9.765625e-04  1.4142135623734411  2.446204e-13
 8: 4.882812e-04  1.4142135623729988  6.814201e-14
 9: 2.441406e-04  1.4142135623727026  2.775923e-13

なんと、最初の時点でもうほとんど精度の限界に迫った値が出ています。これに続けてさらに order を計算すると、

コードと実行結果

 0: -nan
 1: -nan
 2: -nan
 3: -1.033947
 4: -0.527932
 5: -nan
 6: -nan

はい、終わり!

改めて a[0], a[1], a[2], a[3] の誤差だけを横に並べてみます。

1.038416e-02  8.122897e-06  7.560223e-10  1.020560e-14
2.602133e-03  5.083898e-07  1.182280e-11  2.669157e-15
6.509145e-04  3.178545e-08  1.819737e-13  9.734573e-15
1.627525e-04  1.986761e-09  1.240373e-14  1.334579e-14
4.068961e-05  1.241843e-10  1.287476e-14  3.454203e-15
1.017250e-05  7.749505e-12  3.454203e-15  1.679999e-14
2.543130e-06  4.810763e-13  1.648597e-14  4.600371e-14
6.357828e-07  4.553268e-14  4.553268e-14  2.446204e-13
1.589457e-07  4.553268e-14  2.400671e-13  6.814201e-14
3.973646e-08  2.223251e-13  6.327473e-14  2.775923e-13

最初の5行を見れば加速するほど真の値に近いことが分かります。最も精度が高い a[3][1] は小数第15位までの16桁が真の値に一致しています。これは64ビット浮動小数点数ではほとんど最高の精度の値が出ているということを意味します。

値の依存関係を見ると、a[3][1]a[2][1], a[2][2] に依存し、それらはさらに、a[2][1], a[2][2], a[2][3] に依存し、元の数列 a[0] まで遡ると a[0][1], a[0][2], a[0][3], a[0][4] に依存しています。つまり、a[0] の最初の5個の精度の低い(しかし情報落ちが少ない)値から最高に近い精度の値が計算できたということになります。

どうですか?加速法の凄さが分かって頂けたでしょうか?

まとめ

今回は微分を精度よく計算する方法について考えました。定義式をそのまま計算するよりも中心差分近似を使った法が精度よく計算することを見ました。さらに、リチャードソン補外法によって収束を加速する方法について見ました。加速法を用いると、元の数列の最初の方を見るだけでも高い精度が得られることが分かりました。一方で、桁落ちのために加速した数列を進めることは必ずしもよくないことも分かりました。

リチャードソン補外法を用いて微分の精度を上げる方法はロンバーグ微分公式と呼ばれます。また、リチャードソン補外法を積分の台形公式に適用したロンバーグ積分公式もあります。

今回述べたような加速法は関数の滑らかさなどに依存して得意不得意が決まるので、万能の処方箋ではありません。実際に使うときはそのことに気をつける必要があります。基本的には order が不安定になった時点でそれ以上は望めないと考えるべきでしょう。

今回の話では特別なトリックを使っていないことを示すために敢えてC言語を使いました。値を見ながら実験するためにはC言語は良い選択肢とは言えませんが、GPU などで大規模並列計算をする際は依然として第一選択になると思われます。前者の目的では、汎用スクリプト言語である Python がよく使われていますが、ベクトル・行列演算では歴史と実績のある統計解析ツールの R や多くの数値解析ツールを含む Octave、また、科学技術計算を指向した比較的新しい言語である Julia といったより特化したツールが存在します。もちろん、Matlab が使えるならば言うまでもありません。フリーソフトの中では、筆者に馴染みのある R の方がベクトル・行列計算用電卓としては使い勝手が良いように思われました。FreeMat や Scilab は筆者の環境ではよく落ちたので使用そのものに難がありました。

というわけで、加速法の話しでした。では。

*1:一般的な FPU の実装ではそのようになっていて、C言語の処理系もそれを使うように実装されているのが普通です。

*2:Lewis Fry Richardson : 1881-1953

*3:建部賢弘(たけべ かたひろ/けんこう):1664-1739