(前編はこちら)

(中編はこちら)

計算誤差を減らす方法

では、どうすれば計算誤差を減らし、計算精度を高められるのであろうか?

図11 絶対値の順にソートして加算を行う、階層的に加算を行うなどの方法が考えられる

大きい数にずっと小さい数を足すと誤差が大きくなるので、加算する数値を大きさの順にソートして小さい方から順に足し算を行うという方法が考えられるが、計算に先立ってすべての数が揃っていてソートできるとは限らない。また、X、Y、Zの要素がある重力多体問題のベクトルの加算の場合は、Xの大小でソートしても、Y、Zは大小の順にはならないという問題がある。

単純に、前の結果に次の数を順に足しこむのではなく、ペアとなる2つの数の足し算を行い、次に2つのペアの加算結果を足すというように、ツリー状に足し算をするのは良いアイデアで、合計に至るパスの演算回数が減り、大幅に誤差を減らせるという。

また、順次足しこむ場合の高精度の加算のやり方のいろいろと提案されており、次の図12に示すのは、IEEE 754規格の生みの親とも言われるUCBのKahan教授の方法である。

図12 Kahanの加算。Cという補正項を作り、次の数を加算するときにCを引くことで誤差を減らす

ただし、このアルゴリズムはsumの値が足しこむ数yよりずっと大きいことを仮定しているが、同じ程度の正と負の数がランダムに出てくる場合は、sumはほぼ0で、このアルゴリズムはあまりうまく働かないという。

図13は、Javaがどのハードウェアでも同じように動くようにするために、Gustavsonが提案した計算方法で、AとBの和Cを上位のchと下位のclに分解して計算する方法で、(+)はハードウェアの加算で、○に－はハードウェアの減算を意味する。

右側は有効数字が3桁の10進数を使っての説明で、Aが99.6、Bが1.75のケースで、A+Bを3桁に丸めると上位のchは101となる。A>Bであるので、clは1.75から(99.6－101)=1.4を引き、0.35となる。これは上位が101、下位が0.35で正しい答えになっている。一方、これを逆にやると、下位が0.3となり誤差がでる。

この方法を使って、多くの数の合計を求める場合は、chとclの形で順次合計を求めて行き、最後にchとclの和を計算すればよい。

この方法では、合計をもとめる個数Nが100万以下の場合は、常に同じ答えが得られ、Nが200万以上の場合は多少、バラつきが出る程度で、単純に順次加算していく方法に比べると圧倒的に精度が高いという。ただし、約2倍の計算時間が掛かる。

図13 GustavsonがJava用に提案した方法。AとBの和Cを上位のchと下位のclに分解して計算する

結論としては、誤差はGPUだけの問題ではなく、CPUでも同様に発生する。誤差は近似計算である浮動小数点を使うことで発生する根本的な問題である。

高性能のハードウェアができ、大量の演算ができるようになってきているが、演算回数が増えると誤差も増えるので、プログラマはそれに対処する必要がある。誤差を減らすにはいくつかの対処法がある。また、64ビットの倍精度で計算するのは現実的な解であるという。

図14 結論:誤差はGPUだけの問題ではない。浮動小数点を使うことで発生する根本的な問題。演算回数が増えると誤差も増えるので、プログラマはそれに対処する必要がある。それにはいくつかの対処法がある。また、64ビットの倍精度で計算するのは現実的な解

筆者として少し補足しておくと、図8の10個の数値の最大のものは6.94764E-13で最少のものは-8.50532E-13である。この2つの数値を足すと-1.55768E-13で、元の数の1/4以下となっており、10進では有効桁数は減らないが、2進では2ビット程度有効数字が少なくなってしまう。このように、絶対値が同じ程度の正と負の数を足すと、有効数字が減り、誤差が大きくなってしまう。

この発表の例では32ビットの単精度の浮動小数点演算を使っているので、図1に示したようにfは23ビットであるが、64ビットの倍精度の浮動小数点演算を使えば、fは52ビットとなり、ULPは2の-29乗(約5億分の1)に減少する。従って、倍精度で100ULPの誤差がでたとしても、単精度の場合の1ULPの500万分の1の誤差であり問題にならない。

NVIDIAの科学技術計算用のKepler GPUの場合、倍精度の演算性能は単精度の1/3であり、単精度を使った方が3倍性能が高い。従って、単精度で済む場合に、倍精度で計算するのは無駄であるが、誤差が累積する長い計算を行う場合は倍精度の64ビット浮動小数点演算を使うのが安全である。

また、図13に示したGustavsonの加算アルゴリズムを使うと約2倍の計算が必要となるが、それでも倍精度を使うのと比べて1.5倍の性能となる。また、科学技術計算用ではなくグラフィックス専用のKepler GPUは倍精度専用のハードウェアをもっておらず、倍精度演算の性能は単精度演算の1/24の性能しかない。このため、精度の高い加算を大量に行う場合は、倍精度浮動小数点演算を使うよりもGustavsonの加算アルゴリズムを使う方がよい。

また、この発表では触れられなかったが、それぞれの加算を行うとき、丸めモードとして切り上げを選べば、最大の誤差があってもこの値を超えないという上限値が得られる。また、切り下げを行うと、これ以下にはならないという下限値が得られる。このようにして演算ごとに上限と下限を計算し、次の演算でも正しい上限と下限を計算するようにして行けば、正しい答えは、上限と下限の間に入ることが保証される。このように上限と下限を計算するやり方を区間演算(Interval Arithmetic)という。区間演算を行うと、最終結果の真の値は上限値と下限値の間にあることは保証される。この上限値と下限値の開きが小さい場合は問題ないが、大きく異なっている場合は誤差が大きい計算が含まれていることを示しており、より高精度の演算を行う必要がある。