数値計算の実装と検証 | 数値解析 I

1 検証と妥当性確認

数値実験の誤りには、性質の異なる二つの源がある。両者を混同すると、正しいコードを疑い続けたり、不適切なモデルを精緻化し続けたりする。

定義 1.1. 数理モデル（方程式） $\mathcal{M}$ を離散化・実装して数値解 $A(h)$ を得る過程について、次を区別する。

検証 (verification)（verification）——「方程式を正しく解いているか」。離散解 $A(h)$ が、 $h \to 0$ で $\mathcal{M}$ の厳密解 $A$ に、理論が予言する次数で収束することの確認。数学とコードの問題であり、実世界のデータを必要としない。
妥当性確認 (validation)（validation）——「正しい方程式を解いているか」。モデル $\mathcal{M}$ の解が、実際の観測・実験データと許容範囲で一致することの確認。モデル化の問題であり、外部データを要する。

検証は数値解と厳密解の比較、妥当性確認は数値解（またはモデル解）と現実の比較である。

検証を経ずに妥当性確認へ進むと、離散化誤差とモデル誤差が相殺・累積して見分けがつかなくなる。本記事が主題とするのは、外部データを要さずに実行できる検証の手続きである。

2 収束試験と観測収束次数

検証の中核は、刻み幅を細かくしたときに誤差が理論どおりの速さで減るかを測ることである。

定義 2.1. 基準刻み幅 $h_0$ と細分化比 $r > 1$ （多くは $r = 2$ ）を固定し、刻み幅の列 $h_j = h_0 / r^{\,j}$ （ $j = 0, 1, 2, \dots$ ）に対して離散解 $A(h_j)$ を計算する手続きを収束試験 (convergence study) という。厳密解 $A$ が既知なら誤差 $e(h_j) = A(h_j) - A$ を、未知なら連続する解の差 $A(h_j) - A(h_{j+1})$ を、 $j$ に対して記録する。 $e(h_j)$ を $h_j$ の両対数目盛にとった図（収束プロット (convergence plot)）の傾きが、以下の観測収束次数を可視化する。

離散法の設計時に理論上の収束次数 $p$ （例：局所打切り誤差 $O(h^{p+1})$ の一段法が与える大域誤差 $O(h^p)$ 、中心差分が与える $O(h^2)$ ）が分かっていても、実装のバグや仮定の破れは、この理論次数が観測されないこととして現れる。したがって観測次数を理論次数と突き合わせることが検証になる。

命題 2.2. 誤差が $e(h) = C\,h^p + o(h^p)$ （ $h \to 0$ 、 $C \ne 0$ 、 $p > 0$ ）の形をもつとする。細分化比を $r > 1$ とすると、次が成り立つ。

厳密解が既知で誤差 $e(h)$ が測れるとき、 $p = \lim_{h \to 0} \frac{\log\!\big(e(h)/e(h/r)\big)}{\log r}.$
厳密解が未知でも、三つの格子解 $A(h),\,A(h/r),\,A(h/r^2)$ から $p = \lim_{h \to 0} \frac{1}{\log r}\,\log\!\left(\frac{A(h) - A(h/r)}{A(h/r) - A(h/r^2)}\right).$

証明.(1)を示す。仮定より $e(h) = C h^p + o(h^p)$ 、また $(h/r)^p = r^{-p} h^p$ ゆえ $e(h/r) = C r^{-p} h^p + o(h^p)$ 。 $C \ne 0$ だから $\frac{e(h)}{e(h/r)} = \frac{C h^p\big(1 + o(1)\big)}{C r^{-p} h^p\big(1 + o(1)\big)} = r^p\big(1 + o(1)\big).$ 両辺の対数を $\log r$ で割ると $\log(e(h)/e(h/r))/\log r = p + o(1) \to p$ 。

(2)を示す。厳密解を $A$ とすると $A(h) - A = C h^p + o(h^p)$ 。差をとると $A(h) - A(h/r) = \big(A(h) - A\big) - \big(A(h/r) - A\big) = C\big(1 - r^{-p}\big) h^p + o(h^p),$ $A(h/r) - A(h/r^2) = C\big(1 - r^{-p}\big)(h/r)^p + o(h^p) = C\big(1 - r^{-p}\big) r^{-p} h^p + o(h^p).$ $C(1 - r^{-p}) \ne 0$ （ $r > 1$ より $1 - r^{-p} > 0$ ）だから、両者の比は $\frac{A(h) - A(h/r)}{A(h/r) - A(h/r^2)} = \frac{C(1 - r^{-p}) h^p\big(1+o(1)\big)}{C(1 - r^{-p}) r^{-p} h^p\big(1+o(1)\big)} = r^p\big(1 + o(1)\big),$ 再び対数をとって $p$ を得る。▨

(2)は厳密解を要さないため、実問題の検証で実際に使う形である。 $3$ 格子で推定した $p$ が理論値へ近づかないなら、細分化がまだ漸近域に入っていないか、離散化・実装に欠陥があるかのいずれかである。

3 リチャードソン外挿

誤差が刻み幅のべきで展開できるとき、二つの解を線形結合して主要誤差項を消去できる。これにより高精度解と、その場で使える誤差推定の両方が得られる。

定理 3.1 (リチャードソン外挿). 近似が誤差展開 $A(h) = A + c_p\,h^p + c_q\,h^q + o(h^q) \qquad (0 < p < q,\ c_p \ne 0)$ をもち、主要次数 $p$ が既知とする。細分化比 $r > 1$ に対しリチャードソン外挿値を $A_R(h) := \frac{r^p\,A(h/r) - A(h)}{r^p - 1}$ で定めると、次が成り立つ。

$A_R(h)$ の誤差は一段高い次数 $q$ をもつ： $A_R(h) = A + c_q\,\frac{r^{\,p-q} - 1}{r^p - 1}\,h^q + o(h^q) = A + O(h^q).$
細かい側の解 $A(h/r)$ の誤差は、二つの解の差から主要次数まで推定できる： $A(h/r) - A = \frac{A(h/r) - A(h)}{r^p - 1} + o(h^p).$

証明.(1)を示す。展開に $h \mapsto h/r$ を代入すると $(h/r)^p = r^{-p} h^p$ 、 $(h/r)^q = r^{-q} h^q$ より $A(h/r) = A + c_p\,r^{-p} h^p + c_q\,r^{-q} h^q + o(h^q).$ 両辺を $r^p$ 倍して $A(h)$ を引くと、 $h^p$ の項がちょうど打ち消し合う：

r^p A(h/r) - A(h) = \big(r^p - 1\big) A + \big(c_p h^p - c_p h^p\big) + c_q\big(r^{\,p-q} - 1\big) h^q + o(h^q).

右辺の $h^q$ の係数は $c_q(r^{p-q}-1)$ で、 $r^p - 1 > 0$ で割ると $A_R(h) = A + c_q\,\frac{r^{\,p-q} - 1}{r^p - 1}\,h^q + o(h^q).$ $q > p$ より $r^{p-q} - 1 = r^{-(q-p)} - 1$ は $0$ でない有限値だから、主要誤差は $O(h^q)$ である。

(2)を示す。展開の主要項だけを見ると $A(h) - A = c_p h^p + o(h^p)$ 、 $A(h/r) - A = c_p r^{-p} h^p + o(h^p)$ 。差をとって $A(h) - A(h/r) = c_p\big(1 - r^{-p}\big) h^p + o(h^p) = c_p r^{-p} h^p\,\big(r^p - 1\big) + o(h^p).$ 両辺を $r^p - 1$ で割ると右辺は $c_p r^{-p} h^p + o(h^p) = \big(A(h/r) - A\big) + o(h^p)$ 、すなわち $A(h/r) - A = \frac{A(h) - A(h/r)}{r^p - 1} + o(h^p) = \frac{A(h/r) - A(h)}{-(r^p-1)}\cdot(-1) + o(h^p).$ 符号を整理すれば主張の等式を得る。▨

外挿値 $A_R(h)$ を新たな $A(\cdot)$ とみなして展開の次の項 $c_q h^q$ に同じ操作を繰り返すと、 $O(h^{q})$ 、さらに高次へと段階的に精度を上げられる。積分に対してこれを台形則から系統的に行うのがロンバーグ積分である。

検算（台形則へのリチャードソン外挿）. $\displaystyle\int_0^1 e^x\,dx = e - 1 = 1.7182818285\ldots$ を複合台形則で近似する。台形則の誤差はオイラー・マクローリン展開により偶数べき $T(h) = A + c_2 h^2 + c_4 h^4 + \cdots$ をもつので、 $p = 2$ 、 $q = 4$ 、 $r = 2$ とし、外挿式は $A_R = \big(2^2\,T(h/2) - T(h)\big)/3 = \big(4\,T(h/2) - T(h)\big)/3$ となる。各 $T$ は等比和の閉形式で厳密に評価し、誤差比を二重確認した。

$n$	$h = 1/n$	台形則 $T(h)$	誤差 $\lvert e(h)\rvert$	観測次数 $p$	外挿 $A_R$	外挿誤差
$1$	$1$	$1.8591409142$	$1.409\times10^{-1}$	—	—	—
$2$	$1/2$	$1.7539310925$	$3.565\times10^{-2}$	$1.982$	$1.7188611519$	$5.79\times10^{-4}$
$4$	$1/4$	$1.7272219046$	$8.940\times10^{-3}$	$1.996$	$1.7183188419$	$3.70\times10^{-5}$
$8$	$1/8$	$1.7205185922$	$2.237\times10^{-3}$	$1.999$	$1.7182841547$	$2.33\times10^{-6}$

観測次数（真誤差から $\log(e(h)/e(h/2))/\log 2$ で算出）は $2$ へ収束し、台形則の理論次数 $O(h^2)$ と一致する（命題 2.2）。外挿値の誤差は $5.79\times10^{-4} \to 3.70\times10^{-5} \to 2.33\times10^{-6}$ と、 $1$ 段細分化ごとに約 $16 = 2^4$ 倍で減っており、定理 3.1の予言する $O(h^4)$ を確認する。 $n = 8$ の台形則が $2.24\times10^{-3}$ の誤差にとどまるのに対し、 $n = 4, 8$ からの外挿は $2.33\times10^{-6}$ ——同じ計算量でおよそ $1000$ 倍の精度である。

4 製作解法（MMS）

検証は厳密解 $A$ の存在を前提とするが、実際の問題の厳密解はまず得られない。そこで解を先に決め、それを厳密解にする方程式を逆算する。

定義 4.1. 微分作用素 $\mathcal{L}$ に対する問題 $\mathcal{L}u = f$ を解くコードを検証したいとする。製作解法 (method of manufactured solutions) （method of manufactured solutions, MMS）とは、次の手順をいう。

十分滑らかな関数 $u^\ast$ （製作解 (manufactured solution)）を任意に選ぶ。物理的意味は要らない。
強制項を $f^\ast := \mathcal{L}u^\ast$ と解析的に計算する。これにより $u^\ast$ は $\mathcal{L}u = f^\ast$ の厳密解になる。境界・初期条件も $u^\ast$ から定める。
強制項 $f^\ast$ を入力してコードで数値解 $u_h$ を求め、既知の $u^\ast$ との誤差 $\lVert u_h - u^\ast\rVert$ を格子ごとに測り、観測収束次数を理論次数と照合する。

MMS は $u^\ast$ を自由に選べるため、あらゆる項（拡散・移流・反応・非線形・境界）を励起する解を設計でき、離散化のどの部分が理論次数を出していないかを切り分けられる。例えば区間 $(0,1)$ 上の $-u'' = f$ 、 $u(0) = u(1) = 0$ を検証するには $u^\ast(x) = \sin(\pi x)$ ととれば $f^\ast(x) = -\,(u^\ast)''(x) = \pi^2 \sin(\pi x)$ で、境界条件 $u^\ast(0) = u^\ast(1) = 0$ も自動的に満たされる。これを二階中心差分で解けば、 $L^\infty$ 誤差が $O(h^2)$ で減ることを確認できる。観測次数が $2$ に届かなければ、境界の扱いか行列組み立てにバグがある。MMS はモデルの正しさ（妥当性確認）ではなく、離散化・実装の正しさ（定義 1.1の検証）を試す道具であることに注意する。

5 浮動小数点と再現性

以上の手続きは、同じコードが同じ入力に対して同じ出力を返すことを暗黙に仮定している。有限精度の計算では、この「同じ」が思うほど自明でない。

注意 5.1 (浮動小数点加算の非結合性). 実数の加算は結合的だが、丸めを伴う浮動小数点加算 $\oplus$ （§E20.1 系 3.2）は結合的でない：一般に $(a \oplus b) \oplus c \ne a \oplus (b \oplus c)$ 。IEEE 倍精度で $x = 10^{16}$ 、 $y = -10^{16}$ 、 $z = 1$ をとると、 $10^{16}$ の近傍では表現可能な数の間隔（ULP）が $2$ なので $10^{16} \oplus 1$ は最近接丸めで $10^{16}$ に戻る。したがって $(x \oplus y) \oplus z = 0 \oplus 1 = 1, \qquad x \oplus (y \oplus z) = x \oplus (-10^{16}) = 0,$ と結果が食い違う（後者は $z = 1$ が桁落ちで消える。§E20.3 注意 7.3）。帰結として、総和の評価順序が結果を変える。並列リダクションやベクトル化、スレッド数の変更、コンパイラの再結合最適化は加算順序を変えうるため、同じプログラムでもビット単位で異なる答えを返す。これが数値実験のビット単位再現性を破る主要因である。

順序に依存しない総和が必要なら、補償付き加算（カハンの総和）や、ペアワイズ総和、あるいは整数化した固定小数点での確定的リダクションを使う。いずれも、非結合性という事実（注意 5.1）を迂回するための設計である。

注意 5.2 (再現性の要件). 計算実験が再現可能であるとは、第三者（および将来の自分）が結論を同じ手順で再現できることをいう。次を固定し記録する。

乱数シード——確率的な要素（初期値、モンテカルロ標本、確率的最適化）に用いた擬似乱数生成器と、そのシード。
依存バージョン——言語・ライブラリ・BLAS などの数値カーネルのバージョン。丸めや縮約順序は実装に依存しうる（注意 5.1）。
環境——OS・CPU/GPU・スレッド数・コンパイラと最適化フラグ。並列度が結果に影響しうる。
データと前処理——入力データの版と、前処理・スケーリングの全手順。

再現性には二つの水準を区別する。ビット単位再現性は上のすべてを固定して同一のビット列を得る厳しい要求で、並列環境では縮約順序の固定を要する。科学的再現性は、報告した誤差許容（§E20.2 定義 2.1の前進誤差の範囲）内で結論が再現できればよい、より現実的な水準である。論文・実験ログには、どちらの水準を主張するかを明記する。

6 つまずいたら

検証と妥当性確認を取り違えない. 収束次数が出ない（検証の失敗）のにモデルを疑って項を足したり、逆にコードは正しいのに現実と合わない（妥当性確認の失敗）のに実装を直し続けたりしない（定義 1.1）。検証を通してから妥当性確認へ進む。
観測次数は「漸近域」でしか意味をもたない.命題 2.2は $h \to 0$ の極限の主張である。 $h$ が粗いうちは高次項 $o(h^p)$ が効いて観測次数が理論値から外れる。逆に $h$ を細かくしすぎると、離散化誤差が丸め誤差の床（§E20.1 定理 2.2）を下回り、誤差が減らずに増え始める。収束プロットが折れ曲がる手前の傾きを読む。
リチャードソン外挿は誤差展開を前提とする.定理 3.1は $A(h) = A + c_p h^p + \cdots$ という展開と、正しい $p$ の値を仮定する。展開が偶数べきなのに $p = 1$ と誤ると、消去する項を間違えて外挿がかえって悪化する。台形則の $p = 2$ 、前進オイラーの $p = 1$ のように、理論次数を先に確かめる。滑らかさが不足して展開が成り立たない問題では外挿は効かない。
MMS の強制項は解析的に計算する.定義 4.1で $f^\ast = \mathcal{L}u^\ast$ を数値微分で作ると、そこに離散化誤差が混入し、検証したい離散化と切り分けられなくなる。 $u^\ast$ は手計算（または記号計算）で微分できるものを選ぶ。境界・初期条件も $u^\ast$ に合わせる。
「同じコード・同じ入力なら同じ出力」を仮定しない. 加算の非結合性（注意 5.1）により、スレッド数やライブラリ版が変わるだけで結果はビット単位で変わりうる。再現性を主張するなら、シード・バージョン・環境・データを固定して記録し（注意 5.2）、ビット単位再現性と科学的再現性のどちらを主張するかを明示する。

§E20.44数値計算の実装と検証

1 検証と妥当性確認

2 収束試験と観測収束次数

3 リチャードソン外挿

4 製作解法（MMS）

5 浮動小数点と再現性

前提記事