ComputerVisionまとめの部屋

役に立った情報や調査結果をまとめています

浮動小数点の比較①

浮動小数点の比較に関する投稿 2012年版

ブルース・ドーソン著
ソース:https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
 

 
 
この記事は、私が何年も前に書いた「浮動小数点数の比較」に関する記事をより慎重に検討し、ピアレビューを行ったものである。この記事を読めば、テクニカルな浮動小数点数比較というテーマに関するしっかりしたアドバイスと驚くべき見解を得るだろう。
 
私たちはこのテーマに関して待ち望んでいた結果を得ることができた。ここで、私が浮動小数点演算に関してもつ知識の最重要部分を共有したい。
 
 
 

浮動小数点: 数学は難しい

 
あなたは浮動小数点というものがどんなに気の遠くなるほど大変なものか、想像もできないことだろう。シカゴ発の列車とロサンゼルス発の列車が、いつ衝突するかを計算するのは大変だと分かるかもしれないが、それも浮動小数点演算に比べれば取るに足らないものである。
 
 
冗談はさておき、私は浮動小数点演算に含まれる微妙さや意味合いに悩まされるたび、自分が間違っていたことや、検討し損ねた外部因子があったことを発見する。浮動小数点演算は「常にあなたが思うよりも複雑である」ということは、覚えておくべき教訓である。私たちはこの記事で浮動小数点数の比較をトピックにするが、この投稿には「余り」の部分があることを心に留めておいてほしい。またこの投稿はいくつかの技術的な提案を与えるが、万能な解決策(silver bullets:銀の弾丸)を与えるものではないことを理解しておくとよい。
 
 
 
この記事は長いシリーズの第五章である。以下に示すシリーズの最初の部分は、この章を理解するためにとても重要であるから参考にしてほしい。 
 

  1. 浮動小数点のトリック
  2. 浮動小数点のバカなトリック ー整数表現インクリメント
  3. 浮動小数点数型で保持するな
  4. 確かに等しそうだ ーVisualStudioの浮動小数点欠陥について
  5. 浮動小数点数の比較 2012年版(本記事)
  6. 0から100桁の浮動小数点精度 ー有効桁数について
  7. C++11 std::async FastFloatFormatFinding ーわずか数分の全テスト
  8. 中間浮動小数点精度 ー驚くほど複雑な式の評価方法
  9. 浮動小数点の複雑性 ー浮動小数点数のよくある罠
  10. 例外的な浮動小数点 ー面白い上にタメになる
  11. 非正常な奇数浮動小数点数のパフォーマンス ー無限大NaNと非正規化
  12. doubleとfloatは比較するな
  13. 浮動小数点精度の再訪:9桁浮動小数点数の移植性
  14. 浮動小数点決定論 ービット単位で一致を得るには?
  15. 40億個の浮動小数点数全テスト!
  16. 円の円周を求めよ ーC++, 定数、浮動小数点の交差点
  17. インテルの誤差範囲10e+18は過小評価である!

 
 
 

等式による比較

 
浮動小数点数の演算は、厳密ではない。
 
例えば「0.1」のような単純な値は、単精度浮動小数点型で厳密に表現することができない。浮動小数点数の精度が限られるということは、演算順序や、中間計算値の精度によっても結果が変わってくることを意味する。またこのことは、2つの浮動小数点数が等しいかどうかの比較が、必ずしも期待通りにはならないことを意味する。
 
GCCコンパイラではこれに関して以下のような警告を示している。
 

警告: 浮動小数点数を「==」や「!=」によって比較することは安全ではありません。

 
 
以下に例を示す。
 
 

float f = 0.1f;
float sum;
sum = 0;
for (int i = 0; i < 10; ++i)
    sum += f;
float product = f * 10;
printf("sum = %1.15f, mul = %1.15f, mul2 = %1.15f\n",
        sum, product, f * 10);

 
 
このコードは「1」を3つの異なる方法で計算している。
繰り返し加算と、2つの異なる乗算だ。3つの異なる結果のうち「1.0」になるのはそのうちのひとつだけである。
 
 

sum = 1.000000119209290,
mul = 1.000000000000000,
mul2 = 1.000000014901161

 
 
ただしこの結果はコンパイラコンパイラ設定に依存する。このことは実際に誤差原因の解明に役に立つ。
 
 
 
ここで何が起こっているのかを見ていこう。
どの結果が正しいと言えるだろうか?
そもそも「正しい」とはどういうことだろうか?
 
 
 
話を進める前にまずは、「0.1」と「float(0.1)」と「double(0.1)」の違いを明確にしておかなけらばならない。
 
C/C++の記述では「0.1」と「double(0.1)」は同じものを指すが、以下では "0.1" とテキストで書いたときには厳密な10進数であることを意味するものとする。一方で、float(0.1)またはdouble(0.1)と書いたときは0.1が表現可能な数に丸められた数値であることを意味するものとする。
 
ここで、float(0.1)とdouble(0.1)は同じ値にはならないことを明確にしておく。float(0.1)のバイナリ表現での桁数はdouble(0.1)よりも少なく、すなわち誤差が大きくなる。
訳注:floatは4byte=32bit、doubleは8byte=64bitなので、それぞれ32桁と64桁のバイナリ(0/1)表現で表されている。
 
0.1の厳密な値、float(0.1)、double(0.1)の実際の値は以下のようになる。
 

Number Value
0.1 0.1
float(.1) 0.100000001490116119384765625
double(.1) 0.10000000000000000555111512312578270211815834045 41015625

 
 
この結果を踏まえて、上記のソースコードの内容を見ていこう。
 
 

1. sum = 1.000000119209290:sumの計算では、float型で表現可能な数値に丸められた数から始まって、10回の加算を繰り返している。そのため内部での各演算結果はfloat型に丸められ、その度に誤差が入り込むことになる。最終的に得られる結果は、1.0ではなく、10 * float(0.1)でもない。しかし得られる結果は、1.0に最も近い1.0以上のfloat型で表現可能な数値であり、非常に1.0に近い数になる。
 
 
2. mul = 1.000000000000000:mulの計算では、float(0.1)の丸められた数値に10を掛けたものになる。そのためエラーが入り込む余地は少ない。そして0.1からfloat(0.1)への切り上げが発生するのだが、ここでは先に10を掛けてからfloatに丸めることによって切り捨てが発生していて、この例のように「丸め」が2回発生すると正しい値を返すこともあることが分かる。「誤った原因」により「正しい結果」を得ている。とはいえ実際の値はfloat(0.1)の10倍にはなっていないので、そういう意味では「誤った答え」と言える!
 
 
3. mul2 = 1.000000014901161:この計算は丸められた数値から始まって、倍精度で10を掛けているので、丸めによる誤差が入り込む余地が少ない。そのため異なる「正しい結果」、つまりdouble型で保持された10 * float(0.1)の、より厳密な値が得られたことになる。
 
 
さて、1の答えは間違っているがfloatひとつ分だけ離れた数値が得られた。2の答えは合っているが(しかし厳密ではない)、3の答えは完璧に合っている(しかし正しくはない)。
 
 

それでどうする?

 
ここで間違った答えをいくつか得ることができたが(倍精度の答えについては無視することにする)、それでどうすればよいのだろうか?
 
我々は、1.0に等しい答えを得たいのと同時に、1.0に「十分近い」と言えるような答えを得たいとも言えるのではないだろうか。
 

イプシロン値による比較

  
浮動小数点数を等式によって比較するのは良くない考えだ。その代わりに、誤差がある閾値またはイプシロン値以内であるかどうかでチェックするのはどうだろうか?
 

 bool isEqual = fabs(f1 ? f2) <= epsilon;

 
この式によって2つの浮動小数点数が「十分に近い」という概念を表現することができ、このときに2つの数が等しいとみなすことができる。
 
では、イプシロン値として適当な値はなんだろうか?
 
 
我々のこれまでの実験結果は、イプシロン値として使用すべき値はおよそ1.19e-7fであることを示している。実際にfloat.hの定義にはこの厳密な値がFLT_EPSILONという名前で定義されている。
 
 
 
これは明確なことである。ヘッダファイルの神様はFLT_EPSILONが唯一で真のイプシロン値であることを語っている!
 
 
 
それ以外のことはゴミである。1.0*FLT_EPSILONと2.0*FLT_EPSILONの差は、隣り合う浮動小数点数の差を表わしている。しかし1.0より小さい値については、イプシロン値としてのFLT_EPSILONはとても大きな値となり、十分小さな数では比較しようとしている値よりもFLT_EPSILONのほうが大きくなる。
 
 
2.0よりも大きい数では隣り合う浮動小数点数の差は大きくなり、FLT_EPSILONを用いて比較するのはやり過ぎで、あまり明確でないチェックになる。つまり、2.0よりも大きな2つの異なる浮動小数点数の差はFLT_EPSILONよりも大きくなることが保証されている。
 
また、16777216よりも大きくなると浮動小数点数比較に用いる適切なイプシロンの値は1.0よりも大きくなる。ここではFLT_EPSILONは意味をなさないし、望むこともない。
 
 
 

相対的なイプシロン比較

 
 
相対的なイプシロン比較の考え方は、2つの異なる数の差が、どの程度の「度合い」であるかを知るために考え出された。矛盾のない結果を得るためには、2つの数の差と、大きい方の数との比を常に考える必要がある。文章で書くとこういうことだ。
 
 

浮動小数点f1とf2を比較するために、差の絶対値 diff = fabs(f1-f2) を計算する。
もし、diffがabs(f1)とabs(f2)の大きいほうのn% よりも小さければ、f1とf2を等しいとみなす。

 
 
これをコードで書くとこういうことになる。
 
 

bool AlmostEqualRelative(float A, float B,
                         float maxRelDiff = FLT_EPSILON)
{
    // Calculate the difference.
    float diff = fabs(A - B);
    A = fabs(A);
    B = fabs(B);

    // Find the largest
    float largest = (B > A) ? B : A; 

    if (diff <= largest * maxRelDiff)
        return true;

    return false;
}

 
 
この関数は悪くない。ほとんどの場合に機能する。制約については後で述べるが、まずはこの記事のポイントとして、私が何年か前に提案したテクニックについて述べたい。
 
 
 
浮動小数点の相対的比較を行うとき、maxRelDiffを、FLT_EPSILONまたはFLT_EPSILONを数倍した値にすると非常によく機能する。それよりも小さな値にすると、イプシロン値にはならない危険性がある。もしもっと大きな誤差が期待されるのであれば、確かに大きなイプシロン値を使用することもできるが、正気を失ってはいけない。しかしながら、適切な「最大相対誤差値」を選択することは、少し複雑で面倒な、かつ非自明なことであり、悲しいことに浮動小数点数型と直接的な関係はない。
 
 
 

ULP、彼は神経質に言った

 
 
隣り合う浮動小数点数隣接する整数表現を持っていることを以前に述べた。このことはつまり、2つの異なる整数表現数の差を見れば、浮動小数点空間で2つの数がどれだけ離れているかを知ることができることを意味する。このことから以下の結論が導かれる。
 
 
 

ドーソンの自明な結果論:
 
同じ符号の2つの浮動小数点の整数表現があるとき、差の絶対値は「2つの数の間にある表現可能な浮動小数点数に1を足したもの」になる。

 
 
 
これを言い換えれば、「2つの整数表現の数の差が1であれば、2つの浮動小数点数は、等しくなく、かつ最も近い浮動小数点数である」ことが言える。差が2であれば、2つの浮動小数点数はやはり近いが、間にひとつの浮動小数点数がある」ことが言える。2つの整数表現の差を見ることによって、それらの間にどれだけの単位の浮動小数点数があるかどうかを知ることができる。
 
このことは「2つの数は2つのULP(Unit in the Last Place)だけ違う」というような表現でまとめることができる。