位置: 数学編 第31回 / 前提: M2(テイラー展開) / 数学編: M31 / この回でできるようになること: 差分近似の誤差の次数を評価し、刻み幅の適切な目安を出すことができる
目次
0. 物語 ― 刻みを小さくしすぎると誤差が増える
デリバティブのリスク管理チームで、オプション価格の感応度を計算する仕事があります。価格モデルから直接微分を取れないので、入力を少しだけずらして差分を取る方法が使われます。「ずらす幅を小さくすればするほど正確になるはずだ」と考えて、刻み幅を極端に小さくしたレポートが上がってきたとします。
ところが結果は期待と正反対です。たとえば自然対数の底の値、およそ2.718を真の答えとして、刻み幅を一億分の一にしたときの誤差はおよそ1.5億分の1でした。ここまでは刻みを細かくするほど誤差が減っています。しかし刻み幅をさらに十億分の一、一兆分の一と縮めていくと、誤差は逆に膨らみ始めます。一千兆分の一まで縮めると、誤差は0.39にまで跳ね上がります。真の値が2.718ですから、答えの14パーセントもずれてしまう計算です。
原因は計算機の内部にあります。コンピュータは数値をおよそ16桁の精度で保持しています。二つの近い数を引き算すると、上の桁が打ち消し合って、残った下の桁には丸めの影響しか残りません。刻みを小さくするほど引き算の結果が小さくなり、その小さな数を刻み幅で割るので、丸めの影響が増幅されるのです。
つまり、刻み幅には「ちょうどよい大きさ」があります。打ち切りの誤差と丸めの誤差がつり合う点が最適で、それより細かくしても粗くしても精度は落ちます。この回では、その最適な刻み幅がどこにあるのかを、式で求めます。
1. 動機 ― 微分を数値で近似するといくら外れるか
この回を読み終えると、差分近似の誤差の次数を評価し、刻み幅の適切な目安を出すことができます。
導関数の値が欲しいとき、解析的に微分できれば話は終わりです。しかし金融の現場では、オプション価格をモデルの出力として数値的にしか得られない場面が多く、微分も差分で近似するしかありません。
差分近似には複数の方法があり、精度が違います。刻み幅 h を小さくしたとき、誤差が h に比例して減るのか、h2 に比例して減るのかは実用上の大きな差です。h を 1/10 にしたとき、前者は誤差が 1/10 になるだけですが、後者は 1/100 に縮まります。
さらに、「h を小さくすればするほど正確になる」という直感は、計算機の世界では途中で裏切られます。浮動小数点の丸め誤差が効いてくるからです。今回は3つの差分公式を定義し、テイラー展開で誤差の次数を証明し、最適な刻み幅の存在まで示します。
この記事で証明すること
- THM-M31-1 前進差分と後退差分の誤差は、刻み幅の一次 O(h) に比例します。
- THM-M31-2 中心差分の誤差は、刻み幅の二次 O(h2) に比例します。二階導関数の標準的な差分近似も同様です。
- THM-M31-3 浮動小数点で計算すると、刻み幅を小さくしすぎたとき丸め誤差が 1/h に比例して増え、全体の誤差を最小にする刻み幅が存在します。
2. 設定と仮定 ― 刻み幅と3つの差分
f は必要な回数だけ連続微分可能な関数とします。点 x での導関数を差分で近似するために、刻み幅 h>0 を決めます。
- 仮定A1: f は x の近傍で、各定理が要求する回数だけ連続微分可能(最大で C5 級、五階まで連続微分可能)。
- 仮定A2: 計算は倍精度浮動小数点(IEEE 754 binary64)で行う。計算機イプシロンは ϵm≈2.22×10−16。
- 仮定A3: 数値例では h を 10 の累乗で変える。
3つの差分公式を定義します。
前進差分:
D+f(x)=hf(x+h)−f(x)(2.1)
後退差分:
D−f(x)=hf(x)−f(x−h)(2.2)
中心差分:
D0f(x)=2hf(x+h)−f(x−h)(2.3)
さらに、二階導関数の標準的な差分近似も定義します。
D2f(x)=h2f(x+h)−2f(x)+f(x−h)(2.4)
以上が、本記事で扱う近似公式です。
3. 主張 ― 中心差分は一段精度が高い
定理 THM-M31-1(前進差分と後退差分の誤差)
f が x の近傍で C3 級であるとき、
D+f(x)−f′(x)D−f(x)−f′(x)=O(h),=O(h)である。具体的には、前進差分の誤差は 2hf′′(x)+O(h2) である。
定理 THM-M31-2(中心差分の誤差)
f が x の近傍で C5 級であるとき、
である。具体的には、誤差は 6h2f′′′(x)+O(h4) である。
また、f が C4 級であるとき、二階導関数の差分近似 (2.4) についても
である。
定理 THM-M31-3(丸め誤差との兼ね合い)
浮動小数点演算で前進差分を計算するとき、f(x) の評価に大きさ ϵm∣f(x)∣ 程度の丸め誤差が乗るとすると、全体の誤差は
Etotal(h)≈2h∣f′′(x)∣+h2ϵm∣f(x)∣であり、これを最小にする刻み幅
が存在する。中心差分についても同様に最適な刻み幅が存在する。
3つの定理を合わせると、差分近似には「刻み幅を小さくすると打ち切り誤差は減るが、丸め誤差は増える」というトレードオフがあることがわかります。
証明を読む(読まなくても記事は完結します)
4. 証明 ― テイラー展開で誤差項を比べる
この節は読まなくても記事は完結します。
この証明を追いたい人は、数学編 第2回「テイラー展開」へ。
THM-M31-1 の証明
テイラー展開(M2)を使います。f が C3 級であるとき、剰余を O(h3) で押さえられるので、
f(x+h)=f(x)+hf′(x)+2h2f′′(x)+O(h3)(4.1)
です。これを前進差分 (2.1) に代入します。分子に現れる f(x) を相殺するため、f(x+h)−f(x) を計算します。
f(x+h)−f(x)=hf′(x)+2h2f′′(x)+O(h3)
両辺を h で割ると、差分の定義に一致する形になります。
D+f(x)=hf(x+h)−f(x)=f′(x)+2hf′′(x)+O(h2)
したがって
D+f(x)−f′(x)=2hf′′(x)+O(h2)
であり、誤差の主要項は h の一次に比例します。
後退差分についても同様に、f(x−h) をテイラー展開します。
f(x−h)=f(x)−hf′(x)+2h2f′′(x)+O(h3)(4.2)
f(x)−f(x−h) を計算すると
f(x)−f(x−h)=hf′(x)−2h2f′′(x)+O(h3)
h で割って
D−f(x)=f′(x)−2hf′′(x)+O(h2)
誤差の主要項は −2hf′′(x) で、符号が前進差分と逆ですが、大きさは同じ O(h) です。 ∎
THM-M31-2 の証明
中心差分の誤差が一段高い精度を持つことを示します。f が C5 級であるとき、テイラー展開をもう一項先まで書きます。h4 の項は f(x+h) と f(x−h) で同じ符号を持ち、剰余は O(h5) です。
f(x+h)=f(x)+hf′(x)+2h2f′′(x)+6h3f′′′(x)+O(h4)(4.3)
f(x−h)=f(x)−hf′(x)+2h2f′′(x)−6h3f′′′(x)+O(h4)(4.4)
中心差分の分子 f(x+h)−f(x−h) を計算します。(4.3) から (4.4) を引くと、偶数次の項が消えます。これが中心差分の精度が高い理由です。
f(x+h)−f(x−h)=2hf′(x)+3h3f′′′(x)+O(h5)
f(x) の項と 2h2f′′(x) の項が打ち消し合っています。2h で割ると
D0f(x)=f′(x)+6h2f′′′(x)+O(h4)
したがって
D0f(x)−f′(x)=6h2f′′′(x)+O(h4)(4.5)
であり、誤差は O(h2) です。前進差分と後退差分の一次の誤差が、引き算によって打ち消されたのです。
二階導関数の差分近似についても示します。(4.3) と (4.4) を今度は足します。奇数次の項が消えます。
f(x+h)+f(x−h)=2f(x)+h2f′′(x)+O(h4)
移項して h2 で割ると
D2f(x)=h2f(x+h)−2f(x)+f(x−h)=f′′(x)+O(h2)
こちらも O(h2) の精度です。 ∎
THM-M31-3 の証明
実際の計算機では、f(x) の値を正確に計算できるとは限りません。浮動小数点で f(x) を評価したときの値を f^(x) と書くと、丸め誤差は
程度です。ここで ϵm は計算機イプシロン(倍精度で約 2.22×10−16)です。
前進差分を浮動小数点で計算したときの全体の誤差を評価します。打ち切り誤差(テイラー展開の切り捨て)と丸め誤差の2つが加わります。
打ち切り誤差は THM-M31-1 の結果から 2h∣f′′(x)∣ 程度です。
丸め誤差は、分子 f^(x+h)−f^(x) の誤差が最大 2ϵm∣f(x)∣ 程度(2回の関数評価それぞれに ϵm∣f∣ 程度の誤差が乗る)であり、これを h で割るので
程度の誤差になります。2つの誤差を合わせると
Etotal(h)≈2h∣f′′(x)∣+h2ϵm∣f(x)∣(4.6)
です。第一項は h に比例して減少し、第二項は h に反比例して増加します。h で微分してゼロとおきます。
dhdE=21∣f′′(x)∣−h22ϵm∣f(x)∣=0
h2 について解くと
h∗=2∣f′′(x)∣ϵm∣f(x)∣(4.7)
が最適な刻み幅です。∣f(x)∣≈∣f′′(x)∣ の場合は h∗≈2ϵm≈3×10−8 程度になります。
中心差分の場合も同様に考えます。打ち切り誤差は 6h2∣f′′′(x)∣、丸め誤差は hϵm∣f(x)∣ なので
Etotal(h)≈6h2∣f′′′(x)∣+hϵm∣f(x)∣
h で微分してゼロとおくと h∗=(∣f′′′(x)∣3ϵm∣f(x)∣)1/3 が得られます。∣f∣≈∣f′′′∣ の場合は h∗≈(3ϵm)1/3≈10−5 程度です。 ∎
5. 反例 ― 刻みを小さくしすぎると悪化する
CEX-M31-1: 丸め誤差による悪化
f(x)=ex の x=1 での導関数を前進差分で近似するとき、h を 10−1 から 10−15 まで小さくしていきます。理論上は h を小さくするほど誤差は減るはずですが、実際にはそうなりません。
h=10−8 で前進差分の誤差は最小の 6.60×10−9 に達し、それより小さい h では丸め誤差が支配的になって誤差が増え始めます。h=10−15 では誤差が 3.90×10−1 にまで悪化します。これは真の値 e=2.7182818285… に対して約 14% もの相対誤差です。
中心差分でも同様で、h=10−5 付近で最小の 5.86×10−11 を記録し、その後は増加に転じます。
このように、刻み幅を小さくすればするほど正確になるという理解は、浮動小数点演算のもとでは正しくありません。
CEX-M31-2: 滑らかでない関数
f(x)=∣x∣ の x=0 での導関数を考えます。∣x∣ は x=0 で微分不可能(左微分が −1、右微分が 1)です。
前進差分は任意の h>0 に対して
D+∣x∣x=0=h∣0+h∣−∣0∣=1
と右微分を返します。後退差分は −1(左微分)を返します。
中心差分は
D0∣x∣x=0=2h∣h∣−∣−h∣=2hh−h=0
となり、h をどれだけ小さくしても 0 のままです。左右の微分の平均を返していますが、この値はどちらの微分とも一致しません。O(h2) の精度は、テイラー展開が成り立つ(つまり関数が十分滑らかだという)ことが前提です。その仮定が崩れると、中心差分の優位性は失われます。
CEX-M31-3: 境界では中心差分が使えない
定義域の端点では f(x+h) または f(x−h) が定義域の外に出るため、中心差分を構成できません。例えば f が [0,1] で定義されているとき、x=0 では f(−h) が存在しません。この場合は前進差分を、x=1 では後退差分を使うしかありません。
中心差分は常に前進差分より良い、という理解は誤りです。境界条件が課された問題では、片側差分が必要な場面が必ず出てきます。
6. 数値例 ― 刻みを変えて誤差の傾きを測る
f(x)=ex の x=1 での導関数の真の値は f′(1)=e=2.7182818285… です。刻み幅を変えたときの誤差を並べます。
| h |
前進差分の誤差 |
中心差分の誤差 |
| 10−1 |
1.41×10−1 |
4.53×10−3 |
| 10−2 |
1.36×10−2 |
4.53×10−5 |
| 10−3 |
1.36×10−3 |
4.53×10−7 |
| 10−4 |
1.36×10−4 |
4.53×10−9 |
| 10−5 |
1.36×10−5 |
5.86×10−11 |
| 10−6 |
1.36×10−6 |
1.63×10−10 |
| 10−7 |
1.40×10−7 |
5.86×10−11 |
| 10−8 |
6.60×10−9 |
6.60×10−9 |
| 10−10 |
1.55×10−6 |
6.73×10−7 |
| 10−15 |
3.90×10−1 |
1.68×10−1 |
両対数での傾き: h=10−1 から 10−4 の範囲で、前進差分の誤差を両対数グラフに載せたときの傾きは平均 1.0049 です。h を 1/10 にすると誤差が 1/10 になっており、O(h) を確認できます。同じ範囲で中心差分の傾きは平均 2.0001 です。h を 1/10 にすると誤差がほぼ 1/100 に減っており、O(h2) を確認できます。
最適な刻み幅: 前進差分の誤差が最小になるのは h=10−8 付近で、誤差は 6.60×10−9 です。f(1)=f′′(1)=e なので (4.7) の理論値は h∗=2ϵm≈2.98×10−8 で、同じ桁に収まっています。中心差分では h=10−5 付近が最適で、誤差は 5.86×10−11 です。理論値は h∗=(3ϵm)1/3≈8.73×10−6 です。
∣x∣ の原点: 前進差分は常に 1、後退差分は常に −1、中心差分は常に 0 です。h をどれだけ小さくしても変わりません。
数値実験は理論と完全に整合しています。
7. 次回へ ― 乱数生成
今回は差分近似の誤差を定量的に評価しました。O(h) と O(h2) の違いは、本編第98回で偏微分方程式を数値的に解く際に直接効いてきます。
次回 M32 では、一様乱数から正規乱数を生成する方法を扱います。モンテカルロ法の基礎になる話題です。
よくある質問
なぜ中心差分のほうが正確なのですか
両側から挟むことで、一次の誤差が打ち消し合うからです。テイラー展開の偶数次の項が引き算で消えるため、残る誤差が二次の項から始まります。
刻みはどこまで小さくすればいいですか
小さすぎると計算機の丸め誤差が効いてきます。前進差分なら h≈ϵm≈10−8、中心差分なら h≈ϵm1/3≈10−5 が目安です。
O(h2) とは何ですか
刻みを半分にすると誤差がおよそ4分の1になる、という意味です。正確には、h→0 のとき誤差が Ch2 以下に収まる定数 C が存在するということです。
参考文献
- Burden R. L. and Faires J. D.、Numerical Analysis、第9版、Brooks Cole、2010年
- Higham N. J.、Accuracy and Stability of Numerical Algorithms、第2版、SIAM、2002年
数値例は compute.py により倍精度浮動小数点で計算し、assert で検算済み。