2 誤差について
数値計算プログラムを作成するにあたって「誤差」についてまず知っておく必要があります。 ここでは実際にまず誤差を体感していただきます。
2.1 ★ 演習課題:誤差を体感する
2.1.1 ステップ1
以下に示すプログラムの実行結果がどうなるか推測して下さい。program gosa implicit none integer i real delta, time delta = 0.1 time = 0 do i = 1, 1000 time = time + delta end do print '("time = ", f18.13)', time end program gosa
2.1.2 ステップ2
その後このプログラム [gosa.f90](誤差を体感するための演習プログラム例) を実際にコンパイル、実行してその答えを確かめて下さい。出力例: time = 99.9990463256836
0.1を1000回足しているので、正しい答えは100ですが、得られた結果には誤差が生じている事がわかります。
2.1.3 ステップ3
その後、5行目の delta = 0.1 の部分を delta = 1.0 としてプログラムを実行してみて下さい。変更前: delta = 0.1 変更後: delta = 1.0
出力例: time = 1000.0000000000000
正しい答えは1000で、得られている結果も1000なので、今回は誤差が生じていません。
2.1.4 ステップ4
プログラムを倍精度に書き換えて(すべての単精度変数及び単精度定数を倍精度に置き換える)プログラムを実行してみて下さい。 (deltaが0.1の場合と1.0の場合両方で実行を行ってみて下さい。)倍精度に置き換えたプログラムコード例 program gosa implicit none integer i double precision delta, time delta = 0.1d0 time = 0d0 do i = 1, 1000 time = time + delta end do print '("time = ", f18.13)', time end program gosa
delta=0.1d0の場合の出力例: time = 99.9999999999986 delta=1.0d0の場合の出力例: time = 1000.0000000000000
単精度の場合と同様にdeletaが0.1の場合は誤差が見受けられ、1.0の場合は誤差が生じていません。 但しdeltaが0.1の場合に生じる誤差は倍精度の場合の方が小さくなっています。
補足: 倍精度実数の定数表現について
ここで、倍精度実数と単精度実数の定数表現はそれぞれ区別される事に関して補足します。 例えばただ単に 1.23 と記述した場合、これはFortran文法において単精度実数を意味します。 倍精度実数は 1.23d0 と記述する必要があります。
良くある誤りとして、倍精度の変数に単精度の値を代入してしまっている場合があります。 以下の例では倍精度変数 a に単精度実数 1.23 を代入していますが、結果は意図しないもの(1.2300000190734863)となっています。
【誤った記述例 1.23 を使用】 double precision a a = 1.23 print *, a end 出力例: 1.2300000190734863
単精度が求められるところでは単精度、 倍精度が求められるところでは倍精度を指定するよう注意が必要です。
【正しい記述例 1.23d0 を使用】 double precision a a = 1.23d0 print *, a end 出力例: 1.2300000000000000
2.2 演習コード実行結果のまとめ
実行結果例(単精度、delta = 0.1): time = 99.9990463256836 実行結果例(単精度、delta = 1.0): time = 1000.0000000000000 実行結果例(倍精度、delta = 0.1): time = 99.9999999999986 実行結果例(倍精度、delta = 1.0): time = 1000.0000000000000
- deltaを0.1とした場合に誤差が出現したが1.0とした場合には誤差が出現しなかった。
- 誤差が出現した場合において倍精度よりも単精度の場合の方が大きな誤差だった。
2.3 浮動小数点表現と丸め誤差について
コンピュータ内部では浮動小数点表現が利用されています。 浮動小数点表現は以下の図のように、符号部、指数部、及び仮数部から構成されています。 符号部、指数部、仮数部を合わせて32ビット(単精度)や64ビット(倍精度)で数値を表現しています。
それぞれの精度で表現可能な値の範囲と有効桁数は、以下のようになっています。
精度 | 表すことができる範囲 | 有効桁数 |
単精度(32bit) | 10-38~10+38 | 約6桁 |
倍精度(64bit) | 10-308~10+308 | 約15桁 |
このように有限のビット数で無限にある実数を表現することになるため、 正確に表現可能な実数は、上記範囲内であっても、その中のごく一部です。
例えば32ビット単精度実数型において 1.0 前後の値で正確に表現できる数値は以下のとおりです。 (以下の図の緑の点)
単精度実数で表現可能な1.0前後の値(緑の点) |
1.0は単精度の浮動小数点表現で正確に表現できる! |
コンピュータの内部表現ではこの例のように0.999999940395355224609375と1の間及び1と1.00000011920928955078125の間に存在する無数の実数を正確に表すことができません。
正確に表すことができない数値については、 真の値に近い値で且つ正確に表現できる値が用いられます。 そのためどうしても表現上の誤差(丸め誤差と呼ばれる)が発生します。
またコンピュータ内部では10進数ではなく2進数の表現が用いられているため 以下の例のように10進数では仮数部1桁で表現できる数値であっても、内部表現では近似値となる場合もあります。
real A ! 単精度実数を宣言 A = 0.1 ! Aに0.1を代入した(つもり)
実際に変数 A に代入される値は 0.1 ではなく単精度の浮動小数点表現で表すことができる0.1に最も近い値 (すなわち 0.100000001490116119384765625 )が代入されます。
参考図:
単精度実数で表現可能な0.1前後の値(緑の点) |
0.1は単精度の浮動小数点表現で正確に表現できない! |
このような事から、例えば以下の様な思わぬ現象が、簡単に発生してしまうので注意が必要です。
program not_equal implicit none real :: a = 0.3 real :: b = 7 real :: c = 2.1 if( c == a*b ) then print *, 'a x b = c' else print *, 'a x b <> c' end if end program 出力例 a x b <> c[ notequal.f90 ] - イコールにならないプログラム例
2.4 その他の誤差について
以下に丸め誤差以外の代表的な誤差を説明します。2.4.1 情報落ち
コンピュータ内部では有限の桁数で実数が表現されていることは前述したとおりですが、 限られた桁数で絶対値の大きな数と絶対値の小さな数を加減算した際に、 絶対値の小さな数が演算結果に反映されない現象が発生します。これも誤差の一つで情報落ちと呼ばれています。
例えば以下のコード例では100000と0.001を1000回足してその結果を表示していますが、 絶対値の小さい数字0.001の足し算が結果に反映されていません。(正しい答えは100001)
program trailing_loss implicit none integer i real :: a = 100000 real :: b = 0.001 do i = 1, 1000 a = a + b end do print '(F20.5)', a end program trailing_loss 出力例: 100000.00000[ loss-of-trailing-digit.f90 ] - 情報落ちのサンプル
2.4.2 桁落ち
絶対値のほぼ等しい二つの数の差を求めた時に、 その結果の有効桁数が大きく減少してしまう現象(桁落ちと呼ばれる)が発生します。 絶対値が近い値の差の計算過程では、上位の桁がゼロになり演算結果は左詰されます。 この左詰めで空いた(右側の)ビットにはゼロが設定されますが、これにより有効桁数が減少し誤差が生じます。
例えば以下の例では 0.1234 - 0.1231 を単精度で計算させてその結果を出力し、それに引き続く行に 0.0003(単精度の浮動小数表現で 0.0003に最も近い値である 3.000000142492353916168212890625E-04) をそのまま(与えられた書式に合わせて)出力しています。 結果は
0.000300005078 0.000300000014となっていて、桁落ちによる大きな誤差が確認できます。
[ cancellation-error.f90 ] - 桁落ちのサンプル
program cancellation implicit none real a, b a = 0.1234 b = 0.1231 print '(F20.12)', a-b print '(F20.12)', 0.0003 end program cancellation 出力例: 0.000300005078 0.000300000014
2.4.3 打ち切り誤差
打ち切り誤差は計算をある時点で打ち切ってしまったことにより生じる誤差です。 これは数値計算において、与えられた時間内に計算を終える事が求められるためで、 計算を続ければより精度の高い結果が得られるような場合であっても計算を打ち切り、 そこまでに求まった最善の値を結果とすることにより生じる誤差です。
例えばcos(x)をテイラー展開を用いて求める場合、次の無限級数を用いることが可能です。
これを下記のように有限個の項で計算することで打ち切り誤差が生じます。
下記にテイラー展開によりcos(x)を求めるプログラム例を示します。
program taylor implicit none double precision x, mycos, factorial integer i, j print *, "Enter x:" read *, x mycos = 0d0 do i = 0, 3 ! compute(2i)! factorial = 1d0 do j = 2, i*2 factorial = factorial * j end do mycos = mycos +(-1)**i * x**(i*2) / factorial end do print *, mycos end program
2.5 誤差を回避するには?
2.5.1 計算精度を上げる
計算上の誤差は上述のとおり有限ビット数で数値が表現されている事に起因します。 誤差を回避するための最も有効な手段の一つに、計算精度を上げる方法があります。 単精度よりは倍精度、倍精度よりは4倍精度といった具合です。現在Fortranで数値計算プログラムを作成する際に利用可能な精度は主に以下の3つがあります。
Fortran型 | 精度 | 表すことができる範囲 | 有効桁数 |
REAL | 単精度(32bit) | 10-38~10+38 | 約6桁 |
DOUBLE PRECISION | 倍精度(64bit) | 10-308~10+308 | 約15桁 |
REAL(selected_kind(30)) | 4倍精度(128bit) | 10-292~10+308 | 約30桁 |
計算精度を上げる事には主に次の2つのデメリットがあります。 一つはメモリ使用量が増える点でもう一つは演算スピードが下がる点です。
メモリ使用量に関しては、 単精度での使用量を1とした場合、 倍精度ではその2倍、4倍精度では4倍必要となります。
演算スピードに関しては大凡以下の様な事が言えます。 現在利用可能な多く環境では単精度と倍精度の演算はハードウエアにより直接行われ、4倍精度に関してはソフトウエアによるエミュレートで行われます。その結果単精度と倍精度の演算は比較的高速に行うことが可能ですが、4倍精度の演算は単精度や倍精度と比較して一般的に10倍から20倍程度遅くなります。(計算内容やデータサイズ、最適化レベル、コンパイラ等様々な条件に依存するので一概には言えません)
ご参考までに、いくつかの環境で以下のプログラムを実行した場合の比較を示します。
※下記プログラムは3行目のselected_real_kind()の引数を変えることで計算精度を変更できます。(6を与えると単精度、14を与えると倍精度、28を与えると4倍精度となります)
[ precision_perf.f90 ] - 精度の違いによる計算速度を計測するためのサンプルプログラム
program precision_perf implicit none integer, parameter :: wp = selected_real_kind(6) ! 6=SP, 14=DP, 28=QP integer, parameter :: n = 1000 integer t1, t2, t_rate integer i, j, k real(wp), allocatable :: a(:, :) real(wp) val allocate(a(n,n)) a = 1.0_wp val = 1.23_wp call system_clock(t1, t_rate) do i = 1, n do j = 1, n do k = 1, 100 a(j, i) = a(j, i) + val a(j, i) = a(j, i) * val a(j, i) = a(j, i) - val a(j, i) = a(j, i) / val end do end do end do call system_clock(t2) print '(A, F10.3)', 'time it took was:',(t2-t1)/dble(t_rate) end program 出力例: time it took was: 0.958
単精度を1とした場合の実行時間(単位=倍) |
2.5.2 解き方(アルゴリズム)を工夫する
これには、数学的に等価な式で置き換えてみる。計算順序を変えてみる。異なる近似法を採用してみる等が含まれます。ここではあくまで一例としてですが、以下のような一連の20個の数値を足し合わせるプログラムを考えます。
-400000 -800000 -1200000 -1600000 -2000000 -2400000 -2800000 -3200000 -3600000 -4000000 400001 800001 1200001 1600001 2000001 2400001 2800001 3200001 3600001 4000001[ summation_data.txt ] - 入力データ
これらの数値の総和は10となります。
最も一般的な方法では、単純にデータの最初から最後までを順番に足しあわせていく、以下のようなプログラムが考えられます。
program summation implicit none integer i integer, parameter :: n = 20 real total, a(n) ! read numbers open (11, file='summation_data.txt', status='old') read (11, *)(a(i), i=1, n) close (11) ! compute sum total = 0 do i = 1, n total = total + a(i) end do ! print result print '(''true value = '',f5.0)', n/2.0 print '(''computed value = '',f5.0)', total end program[ summation.f90 ] - 単純に足し合わせる例
しかしこのプログラムを実際に走らせて見たところ正しい答えとは程遠い6という結果が得られました。
出力例: true value = 10. computed value = 6.
これは計算過程で、totalの絶対値がある程度大きくなった際に、絶対値の小さな数値が結果に反映されなかった為(情報落ち)です。
一つの改善策としては以下の様に、 totalの絶対値があまり大きくならないように、絶対値の小さい順に数値を並べ替えてから計算させる事が考えられます。
絶対値の小さい順に並べ替えてから… |
program summation_sorted implicit none integer i, j integer, parameter :: n = 20 real total, tmp, a(n) ! read numbers open (11, file='summation_data.txt', status='old') read (11, *)(a(i), i=1, n) close (11) ! sort numbers by magnitude (selection sort) do i = 1, n - 1 do j = i + 1, n if (abs(a(i))>abs(a(j))) then tmp = a(i) a(i) = a(j) a(j) = tmp end if end do end do ! compute sum total = 0 do i = 1, n total = total + a(i) end do ! print result print '(''true value = '',f5.0)', n/2.0 print '(''computed value = '',f5.0)', total end program 出力例(正しい結果が得られた!): true value = 10. computed value = 10.[ summation-sorted.f90 ] - 絶対値の小さい順に並べ替えてから足し合わせる例
計算順序を工夫したことで正しい結果が得られました。
2.5.3 先人の知恵を借りる
複数の数値の総和を求める方法には、実は 並べ替える方法よりも優れていて、より適用範囲の広い方法(アルゴリズム)が知られています。 ここでは詳しく説明をしませんが、誤差を補正しながら数値を足しあわせて行く方法で、 考案者の名を取り「Kahan summation algorithm=カハンの加算アルゴリズム」と呼ばれているものがあります。以下にKahan summation algorithmのコードセグメントを示します。
! Kahan summation algorithm compensation = 0 total = 0 do i = 1, n y = a(i) - compensation t = total + y compensation =(t - total) - y total = t end do
この方法で総和を求めると以下のように並べ替えを行う事なく正しい結果が得られています。
出力例: true value = 10. computed value = 10.[ summation-kahan.f90 ] - Kahan summation algorithmを用いて足し合わせるプログラム例
このように、一般的に良く行われる計算については、既に知られているアルゴリズム(解法)が存在する場合があります。
自分が行いたい計算について、既存のアルゴリズムが存在するかどうか調べる価値は十分あると言えます。
2.6 絶対誤差と相対誤差という考え方
絶対誤差とは近似値と真の値の差です。
絶対誤差=|真の値-近似値| (絶対誤差はただ単に誤差と呼ばれることもあります)
誤差には相対誤差という考え方もあります。 相対誤差は真の値に対する誤差の割合ととらえることもでき、以下の式で求まります。
相対誤差=絶対誤差/|真の値|
ここで一つの例を考えてみます。 ある人の身長の真の値が160cmである場合に測定値として160.1cmが得られたとします。 この場合の絶対誤差は0.1cmです。同様に蟻の場合を考えて見ます。 蟻の身長の真の値が1cmで、測定値として1.1cmが得られた場合も、絶対誤差は同じ0.1cmです。
同じ0.1cmの絶対誤差ですが、蟻の場合の方がより重大な誤差のように思えます。
そこで、それぞれについて相対誤差を計算してみると、 人間の場合は 0.000625(0.1÷160)で、蟻の場合は 0.1(0.1÷1)となり、 蟻の場合の方がかなり大きな値である事がわかります。 相対誤差は真の値の大きさに対する誤差の比率を表す数値なので、このような誤差の重要度を判断する手がかりとして役立てる事が可能です。