数値計算に混入する誤差を紹介する.
注）このほかに現象を数値計算で理解・予測する際にはモデル化誤差という数値計算ではどうにもならない誤差もある.
浮動小数点数を使用した演算に混入する誤差について述べる
浮動小数点数同士の演算（加減乗除など）の結果は，浮動小数点数で表せるとは限らない． 例えば，10進数で仮数部3桁の浮動小数点演算を考え，2/3を計算すると，
$$ 2.00\times 10^0/3.00\times 10^0=0.66666666...\times 10^0 $$となり，仮数部3桁に収まらない．仮数部の4桁目で四捨五入を行うと
$$ 6.67\times10^{-1} $$となる．このときの計算値と真値との差
$$ 6.67\times10^{-1}-6.6666666...\times10^{-1}=3.3333333...\times10^{-4} $$が丸め誤差である．
IEEE 754は2進数浮動小数点数なので，基本的には0捨1入で丸められる． 例えば，10進数の「0.1」をIEEE 754のbinary64に変換してみると
\begin{align*} 0.1_{(10)} &= 0.000110011001100110011..._{(2)}\\ & = 1.10011001100110011...\times 2^{m}, \end{align*}$m$は$-4+1023=1019$の2進表現，つまり$m=01111111011$. 仮数部は無限小数になっているのでそのまま格納出来ない． 小数点以下を 52bit以内とそれ以降で区切って表示すると
$\fbox{1001100110011001100110011001100110011001100110011001}\fbox{10011001100....}$
となり，はみ出た部分の先頭が「1」なので、0捨1入で繰り上げる． 最終的には， 10進数の「0.1」は
$\fbox{0}\fbox{01111111011}\fbox{1001100110011001100110011001100110011001100110011010}$
のように格納されている．
versioninfo()
Julia Version 1.9.0 Commit 8e630552924 (2023-05-07 11:25 UTC) Platform Info: OS: macOS (arm64-apple-darwin22.4.0) CPU: 8 × Apple M2 WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-14.0.6 (ORCJIT, apple-m1) Threads: 2 on 4 virtual cores
x = 0.1;
println(bitstring(x))
0011111110111001100110011001100110011001100110011001100110011010
10進数の0.1は計算機には正確に格納できず，少しだけ0.1より大きい値で格納されており， 丸め誤差が含まれている
$$ \displaystyle\frac{1}{10}\approx{\color{red}{0.1000000000000000}}055511151231257827021181583404541015625 $$using Printf
@printf("%.55f", x)
0.1000000000000000055511151231257827021181583404541015625
さらに, 0.1を逐次的に加えていくと計算結果は丸め誤差の影響により, 不可解な計算結果を返す事がある.
x = 0.0;
for i = 1:100
x += 0.1;
end
@printf("%.55f\n", x)
9.9999999999999804600747665972448885440826416015625000000
IEEE754では，4つの丸めモードを用意している．$\tilde{x}$を実数（$\tilde{x}\in\mathbb{R}$）とする．いま実数 $\tilde x$ が
$$\tilde x=\pm\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\frac{d_2}{2^2}+\dots\right)2^{e}_{(10)}$$であるときに，その近似
$$x=\pm\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\frac{d_2}{2^2}+\cdots+\frac{d_{52}}{2^{52}}\right)2^{e}_{(10)}\in\mathbb{F}$$を採用する丸めは原点方向への丸め（切り捨て）と呼ばれている．
また $\tilde x$ の近似として誤差$|x-\tilde{x}|$の値が最小になる$x\in\mathbb{F}$，すなわち
$$ |x-\tilde{x}|=\min_{y\in\mathbb{F}}|y-\tilde{x}| $$をみたす$x\in\mathbb{F}$を採用する方法を最近点への丸めという．しかし $\tilde x$ が2つの浮動小数点数
$$ x_1=\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\frac{d_2}{2^2}+\cdots+\frac{d_{52}}{2^{52}}\right)2^{e}_{(10)},~~~x_2=\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\frac{d_2}{2^2}+\cdots+\frac{d_{52}+1}{2^{52}}\right)2^{e}_{(10)} $$の中点となると$|x_1-\tilde{x}|=|x_2-\tilde{x}|$が成り立ち，$x$は一つに決まらない．このときは仮数部の最後のビットが0になる方に丸めが実行される．これを最近偶数への丸めという．
eps = 2. ^ (-52) # unit round off
x = 1. + eps/4.
@printf("%.50f\n", x)
println(bitstring(x))
y = 1. + 3*(eps/4.)
@printf("%.50f\n", y)
println(bitstring(y))
z = 1. + (eps/2.)
@printf("%.50f\n", z)
println(bitstring(z))
1.00000000000000000000000000000000000000000000000000 0011111111110000000000000000000000000000000000000000000000000000 1.00000000000000022204460492503130808472633361816406 0011111111110000000000000000000000000000000000000000000000000001 1.00000000000000000000000000000000000000000000000000 0011111111110000000000000000000000000000000000000000000000000000
x = 2.
@printf("%.50f\n", x)
println(bitstring(x))
y = 2. - eps
@printf("%.50f\n", y)
println(bitstring(y))
z = 2. - eps/2.
@printf("%.50f\n", z)
println(bitstring(z))
2.00000000000000000000000000000000000000000000000000 0100000000000000000000000000000000000000000000000000000000000000 1.99999999999999977795539507496869191527366638183594 0011111111111111111111111111111111111111111111111111111111111111 2.00000000000000000000000000000000000000000000000000 0100000000000000000000000000000000000000000000000000000000000000
次は$\tilde{x}$よりも必ず大きな$a\in\mathbb{F}$，すなわち
$$ a=\min\left\{x\in\mathbb{F}:x\ge\tilde{x}\right\} $$をみたす$a\in\mathbb{F}$を採用する方法を$+\infty$方向への丸め（上向き丸め）という．逆に$\tilde{x}$よりも必ず小さな$b\in\mathbb{F}$，すなわち
$$ b=\max\left\{x\in\mathbb{F}:x\le\tilde{x}\right\} $$をみたす$b\in\mathbb{F}$を採用する方法を$-\infty$方向への丸め（下向き丸め）という．これらより
$$ b\le\tilde{x}\le a $$が常に成立する．特に精度保証付き数値計算ではこの2つの丸めモードを利用することで，厳密な包含を得ることができるようになる．まとめると
丸め誤差は実数を浮動小数点数で近似する際の誤差であった．ここでは浮動小数点数を用いた演算のその他の問題点を紹介する．
極めて近い数どうしの減算によって，誤差が著しく大きくなってしまう現象． 2つの浮動小数点数
$$ x=\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\dots+\frac{d_p}{2^p}+\frac{1}{2^{p+1}}+\cdots+\frac{b_{p+2}}{2^{p+2}}+\dots+\frac{b_{52}}{2^{52}}\right)2^{e}_{(10)}, $$$$ y=\left(\frac{1}{2^0}+\frac{d_1}{2^1}+\dots+\frac{d_p}{2^p}+\frac{0}{2^{p+1}}+\cdots+\frac{c_{p+2}}{2^{p+2}}+\dots+\frac{c_{52}}{2^{52}}\right)2^{e}_{(10)} $$が$x>y$とし，仮数部の最初から $p$ ビットが等しいとする．このとき
$$ x-y=\left(\frac{1}{2^0}+\frac{b_{p+2}-c_{p+2}}{2^1}+\dots+\frac{b_{52}-c_{52}}{2^{52-p-1}}\right)2^{e-p-1}_{(10)} $$これよりもともと $52$ 個あった仮数部の情報が，$52-p$ 個に減っている．例を挙げよう．
$b>0$とし，2次方程式 $x^2+bx+c=0$の解の公式
$$ x_1=\frac{-b+\sqrt{b^2-4c}}{2},\quad x_2=\frac{-b-\sqrt{b^2-4c}}{2} $$を考える．いまもしも$b^2\gg c$となるならば，$b$ と $\sqrt{b^2-4c}$が近い数になるので，$x_1$の分子の計算で，桁落ちが起こる．
b = 1e+15;
c = 1e+14;
x1 = (-b + sqrt(b^2 - 4c)) / 2;
println(x1)
println(x1^2 + b*x1 + c)
x2 = 2c / (-b - sqrt(b^2 - 4*c))
println(x2)
println(x2^2 + b*x2 + c)
-0.125 -2.4999999999999984e13 -0.10000000000000002 0.0
絶対値の大きさが極端に違う2 数の加減算を行った時，小さいほうの数値の下位の桁が失われてしまう現象．
println((3.14159265358979+1e10)-1e10)
println(3.14159265358979 + (1e10-1e10))
3.141592025756836 3.14159265358979
println(1e48+543.2-1e48-1e36+123.4+1e36)
println(1e48-1e48-1e36+1e36+543.2+123.4)
println(1e48-1e36-1e48+1e36+543.2+123.4)
0.0 666.6 -3.9230984419161617e30
無限回行うべき計算を有限回の計算で置き換えることにより生じる誤差．計算機は有限回の四則演算しかできない．そのため，無限級数や収束列のような値を求めるためには，有限項で打ち切った近似値を用いる．その際に誤差が生じる．
（例）Taylor展開の打ち切り誤差，Newton法の打ち切り誤差，数値積分の打ち切り誤差
function exp_taylor(x) # naive implementation of "exp" function
s = 0.0;
t = 1.0;
i = 1;
while true
s += t;
if abs(t) < abs(s) * 1e-15
break
end
t *= x / i;
i += 1;
end
return s
end
println(exp_taylor(20.))
println(exp(20.))
println(exp_taylor(-20.))
println(exp(-20.))
4.851651954097903e8 4.851651954097903e8 6.147561828914626e-9 2.061153622438558e-9
離散化誤差の説明はそのうち加える。
本資料も筆者が学生の頃に精度保証付き数値計算を教えて下さった柏木雅英先生の「数値解析特論」の講義資料が基になっています. また, 以下のような文献・Web ページ等を参考にこの文章は書いています.
(Twitterとかでも度々話題に上がる名著. IEEE754 の制定の年にすでに浮動小数点数に対する注意が詰まっている書籍が出版されている. 桁落ち, 情報落ちなどの誤差に詳しい)
(数値解析学の現在最も詳しい教科書)
(精度保証付き数値計算の教科書. 浮動小数点数および区間演算に詳しい. この1章が読めたら大したもの)
(IEEE 754 浮動小数点数を細かく紹介し, 丸め誤差の詳細, および区間演算について触れている)
(数値解析の超有名人によるブログ記事, 丸めについて端的にまとめられている)