誤差が出るケースと簡単な計算例
図4は、丸めにより演算結果に違いが出る原因をまとめたもので、fのビット数が有限であることから、1/3を次々と加えても2.0にはならないPrecisionが問題のケース、有効桁数より小さい値(この例では有効数字は小数点以下2桁で、それに0.00333)を加えても1.0から増加しないBelow ULP(Unit of Least Precision:最下位のビットの重み))のケース、(A+B)+CとA+Bを先に計算した場合とA+(B+C)とB+Cを先に計算した場合のように、計算の順序が違うと答えが異なるケース、丸めの回数で値が異なるDouble Roundingのケースなどがある。また、sin(x)のような関数や1/xのような逆数はIEEE 754規格では定義されておらず、多項式近似で計算するので、ULPで数ビット程度の誤差が出る場合がある。
|
|
図4 演算結果に違いがでる原因 |
1.0の1/1000000を100万回加えると、正しい結果は2.0になるはずである。
|
|
図5 1.0に100万分の1を100万回加算するプログラムとその結果。1.95367…から2.0をわずかに超える値までいろいろな結果が出る |
図5に示すように、1.95367…から2.0をわずかに超える値まで3種類の結果が出る。このうち、どれがGPUの結果であろうか? 実は、これはすべてCPUで計算した結果である。
最初の結果は倍精度の演算を使い、コンパイラで最適化を行わないプログラムで計算した場合の結果である。2番目は-O0で軽い最適化をコンパイラに指示した場合で、最適化なしと同じ答えが得られている。3番目は-O3で最高レベルの最適化を行った場合で、演算の順序が変わり2.0より少し大きい値になっている。
4番目と5番目は倍精度より長い80ビットで演算を行うx87を使って計算を行い、4番目は演算結果を倍精度に丸めた場合で、最後は単精度に丸めた場合の計算結果である。
|
|
図6 最初の結果は倍精度を使い、最適化を行わないで計算した場合、2番目は-O0で軽い最適化を行った場合、3番目は-O3で最高レベルの最適化を行った場合。4番目はx87を使い結果を倍精度に丸めた場合で、最後は単精度に丸めた場合 |
1番目、2番目と4番目は同じ値であるが、正しい答えである2.0より小さい値になっている。3番目の最適化レベルを最高にした場合は、2.0を超える値になっているが、ここでは2.0にもっと近い値となっている。また、最後の単精度に丸めるケースは大きな誤差がでている。
このように、CPUを使っても、単純な加算の場合でもやり方によって、いろいろな答えが出てくる。
重力多体問題の計算
Lars Nyland氏は、重力多体問題の計算を例にとり、浮動小数点演算の誤差がどうなるかと誤差を減らす対策を説明した。多数の星の重力相互作用を計算する場合、ペアの距離と質量の積から重力が計算されるが、これをすべてのペアに対して計算してそれぞれの物体に働く力を計算して合計する方法を直接法という。これに対して、遠方にある物体の重力はまとめて多項式で近似して計算する高速マルチポール法という計算法がある。
|
|
図7 中心の灰色の箱の中の物体に対して、周囲の物体が及ぼす重力を計算する。すべての物体との重力を計算する直接法と、遠方の濃い灰色の箱の中の物体はまとめて扱う高速マルチポール法がある |
右のグラフは物体の数を横軸にとり、縦軸は計算誤差を示しており、ほぼフラットな3組の線は高速マルチポール法の計算で、多項式の項数pが上から順に4、8、12の近似となっている。なお、○の点はCPU、●の点はGPUの計算結果で、FMAを使うGPUの方が、若干誤差が小さい。 そして、右上がりの斜めの線がGPUでの直接法の計算の誤差である。直接法は物体の数が1万以下の場合は、p=12のマルチポール法より誤差が少ないが、物体数が10万を超えると近似計算であるp=12のマルチポール法より誤差が大きくなってしまう。厳密解法が近似計算に負けるというこの結果は筆者にとってはショックであった。
これは直接法では計算量はO(N2)で、中央の箱の中にある1つの点についての計算はNに比例して増加し、誤差が累積するのに対して、高速マルチポール法の計算量はO(N)で1つの物体に働く力の計算量はNによらないので、誤差は物体数によらず一定となるからである。
重力多体問題では、全方向にほぼ一様に物体があり、X、Y、Zそれぞれの＋方向に働く力と－方向に働く力が同じくらい存在するというのが一般的な状況である。この状況を単純化したものが図8で、白抜きの中に書かれた10個の数値を足すという計算を行ってみる。右側の棒グラフは、この10個の数値を順に並べたものである。 図8に示す10個の数をランダムに順序を変えて単精度で加算するという計算を10万回行い、得られた結果のヒストグラムが図9である。
|
|
|
図8 白枠の中の10個の数を単精度で加算すると |
図9 順序を変えて加算すると、47種の異なる答えが得られた。そして、誤差が最も大きいケースでは最下位ビットの25倍というケースがみられる |
図9の横軸は、誤差がULP(fの最下位ビットの重み)の何倍であるかを示している。この図を見ると、+1ulpのケースが最も多いが、最悪のケースは+25ulpとなっており、マイナス側も-23ulpという誤差のケースがある。
図9では10個の数字を加算していたが、これをより一般化して、1000個、2000個、3000個、4000個、5000個とした場合の誤差の分布を示したのが図10である。
|
|
図10 数字の数を1000、2000、3000、4000、5000とした場合の誤差の分布 |
加算する数が多くなるにつれて分布がなだらかになり、正しい答えが得られる確率が低くなる。特に、N=5000になると、50ulpの誤差となるケースの出現頻度も無視できない率になっている。
(後編に続く)
