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)だけ違う」というような表現でまとめることができる。
 
 

人工知能と産業・社会

〜第4次産業革命をどう勝ち抜くか〜

まとめ
これまでの人工知能ブームと、近年話題になっている人工知能はどう違うのか。その応用性や社会に与える影響を、政治家の視点でわかりやすく説明、予測を行っている。人工知能の背景や将来を知る入門書として良い本。

  • 第1次AIブーム:機械学習パーセプトロン、探索木
  • 第2次AIブーム:知識ベースと推論機構(ベイジアンネットワーク)を組み合わせて、問題に対する答えを導き出す。Siri、機械翻訳などのエキスパートシステム
  • 第3次AIブーム:ディープラーニング
  • 探索木のアルゴリズムで解決できるのはごく一部の単純かされた問題のみ。
  • 第2次AIブームにおける課題は、組み合わせが多層になる場合に正解を見つけるための探索の組み合わせが爆発してしまう。効率的に正解を見つける方法が必要だった。

  • 知能とは何か?
    知能とは環境適応能力のことであり、判断力・理解力のこと。学習能力、思考能力はあいまいな定義であり、人工知能の定義もあいまいである。
  • 知能は試行錯誤により構成されるもので、外界との相互作用のなかで洗練されていくもの。適応は動的であり、知識は静的である。

  • 機械学習は、教師あり学習と教師なし学習に分類される。教師あり学習とは、入力に対する正しい出力(訓練データ)が与えられてそれを再現するパラメータを決定する方法。教師なし学習とは、入力データのみが与えられて類似性の近いものをグループ分けする方法。正解はなく、最終的にいくつのまとまりに分けられるかが重要。
  • 強化学習とは、与えられた入力データにたいして出力を生成し、その出力にたいして報酬を与えて、報酬が高くなるような出力を試行錯誤して探す方法。正しい出力を与えられている分けではないので、教師あり学習とは異なる。
  • これまでの人工知能は「概念」(対象のもつ特徴の総和)を理解することはできなかった。すなわちチューリップの花が赤いことを覚えられても、チューリップとはどういうものかを理解することはできないので、様々な花がある中で赤いチューリップを見つけることはできなかった。ディープラーニングによりそれが可能となった。

  • 第3次AIブームを引き起こしたディープラーニングによって、AIが人間の知能を超える「シンギュラリティ(技術的特異点を引き起こす可能性が出てきた。
  • ディープラーニングは教師なし学習に分類される。従来扱われてきたよりも多層のニューラルネットワークを活用し、入力と出力を同じものにする「教師なし事前学習」により各層を学習させいく。

  • ディープラーニングを学習教材に応用した、アダプティブ・ラーニングの導入の社会的メリットは大きい。生徒にとっては、取りこぼしをなくし空いた時間で創造的発想やディベート力を養うことができる。先生は、平均的な生徒に向けた画一的な授業ではなく、個々の生徒の進み具合に合わせた授業ができる(最適化カスタマイズ学習)。またこの仕組みを提供する企業は、能力をもつ人材を効率良く発掘することができるメリットがある。 
人工知能と産業・社会

人工知能と産業・社会

 

 

特徴点を用いた物体認識

OpenCVでは、特徴点抽出・特徴記述・特徴点マッチングについてさまざまなアルゴリズムが実装されているが、共通のインターフェイスが用意されている(OpenCV3.0)。

共通インターフェイスを使えば、異なるアルゴリズムであっても同じ書き方でアプリケーションを実装することができる。

特徴点抽出は以下の3つのインターフェイスで構成される。

DescriptorExtractorインターフェイスを実行する際は、DescriptorExtractor::computeを呼び出す。第1引数に画像、第2引数に特徴点の配列、第3引数に画像の特徴情報を格納するための行列を指定する。

特徴点抽出からマッチングまで

1.画像内のコーナー(キーポイント)を求める
2.キーポイントを特徴付ける特徴量を計算する
3.キーポイント間の類似度を求める
3.類似度が指定した閾値以下であれば対応点とみなす(=マッチング)

特徴量を用いた物体認識

局所特徴量は、スケールの変化・平行移動・回転・隠蔽(オクルージョン)に対して頑強であるべきである。特徴点を使ったマッチングは、多少間違えがあっても、マッチング線が「特定の物体」に集中していれば(得票数が多ければ)きちんと物体が認識できることを利用して、テンプレートマッチングなどの方法と比較して高速化できる利点もある。

クエリとして与えた画像にどんな物体が写っているかを認識するには、物体モデルデータベースを検索して認識する。それが何かを認識するためには、その物体がモデルデータベースにあらかじめ登録されていないと認識できない。(登録されていない物体は何であるかを認識できない)

物体モデルデータベースは、モデルとする物体の画像から特徴量を抽出して物体モデルデータベースに格納し、各画像から抽出してキーポイントの特徴量を物体IDとひもづけて保存する必要がある。

 

特徴点抽出法の種類

Eigen

最小固有値が十分大きいコーナーを検出する方法

Harris

強いコーナーを検出する方法で、回転普遍性はあるがスケール変化普遍性がない

Fast

より多くのコーナーを検出する方法(OpenCV2.0から追加)

SIFT

特徴量を求める適切な範囲を計算してスケール・回転普遍性をもたせた(特許制約あり)

SURF

スケール変化・回転・輝度変化に対して耐性がある(特許制約あり)

STAR

SIFTやSURFと比較して、ローカライズ(エッジ上などノイズの影響を受けやすい部分を除去する処理)が入っている(OpenCV2.2から追加された)

MSER

領域分割の方法を利用。MSCRは改良版。

ORB

回転普遍性がありSURFの10倍、SIFTの100倍と高速である特許制約がなくライセンス的にも使いやすい特徴量。回転には環境だがスケール変化には比較的弱い。

AKAZE

OpenCV3.0で新たに追加された。SIFTやSURFの欠点を解決し、ロバスト性の向上と高速化を図ったもの。特許制約なし。SURFよりも処理速度が遅いが、これはOpenMPを用いた並列化処理が現時点ではなされていないためである。AKAZEは並列化に向いているのでOpenMPの利用により大幅なスピードアップが見込まれる

特徴量の抽出は手法によって求められる特徴量がまったく違うことがあるので、アルゴリズムの中身を理解して用途に応じて適切な手法を選択する必要がある。

NEON

ARMの標準SIMD命令セットAdvanced SIMD(single instruction multiple data:1命令で複数のデータの演算を行うコンピュータの並列化の形態のこと。パック演算、ベクトル演算とも表現する)の通称。

ARMレジスタの演算処理は1命令で1演算の逐次的処理を行うが、NEONレジスタの演算処理は指定したデータサイズにて1命令で複数の演算処理を行う。ハードを利用して高速化できない部分をさらに高速化するためにNEON化する。

 

NEONコプロセッサを使用する場合、以下の3つの方法から選択することができる。

アセンブリ命令は、ARMプロセッサのパイプラインを考慮したプログラミングが必要となり、コーディングが難しいため、自動ベクトル化とNEON組込み関数の使用する方法が推奨される。

ORB(Oriented FAST and Rotated BRIEF)

ORBによる特徴量検出方法は、BRIEFと同じくバイナリコードを用いた高速な特徴点検出が可能であり、BRIEFにはない回転不変性を導入したものである。

ORBは、BRIEFと同じく画像を多段階に縮小したピラミッド画像に対してFASTによる特徴点検出方法を用いており、スケール不変性をもつ。ただし、オリエンテーションの算出がBRIEFとは異なる。(BRIEFでは規則的にならんだ画素値の長距離ペアの勾配からオリエンテーションを求めていた)

ORBではオリエンテーション方向を、特徴点を中心としたパッチ内の輝度のモーメントを算出し、パッチの輝度重心の方向ベクトルとして求めている。この方法を「Oriented FAST」呼ぶ。求めたオリエンテーションにしたがって輝度値を規格化することにより、ORBは回転不変性を得ている。

ORBでは、パッチ内の観測点をランダムに決めるBRIEFと同じ方法を採用するが、さらに高精度なマッチングを行うために、学習を用いてパッチ内の観測点を決定している。観測点決定の基準は以下のようになっている。

(1)バイナリコードのビット分散が大きいこと 
⇔ ある観測点のペアに注目したとき、特徴点ごとにバイナリコードができるだけばらつくようにピクセルを選ぶこと 

(2)バイナリコードのビット間の相関が小さいこと
⇔ ある観測点のペアと別のペアを比較したときに、バイナリコードができるだけ一致しないペアを残すこと

画像がオリエンテーション方向に規格化されるので、上の基準で選択された観測点の並びを見ると、(1)のビット分散の大きいペアだけを選択すると、縦に大きく並んだペアが高い頻度で存在するようになっている。また、(2)のビット間相関が小さいものも含めてペアを選択すると、縦に並んでいるペアがだいたい等間隔になっている、という結果が見られる。

 

https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0CCAQFjAAahUKEwiyo8mc-77IAhXj46YKHUWDCzM&url=https%3A%2F%2Fwww.willowgarage.com%2Fsites%2Fdefault%2Ffiles%2Forb_final.pdf&usg=AFQjCNHRzIFBK9unVZdtJPg7WpSXkbucKA

 

BRISK(Binary Robust Invariant Scalable Keypoints)

以前の日記で、スケール・回転不変な特徴点検出方法としてSIFTやSURFについて述べた。SIFTの課題を解決したのがSURFであり、DoG画像生成や勾配ヒストグラム生成の計算コストが高い問題を、SURFでは積分画像を利用することにより10倍の高速化を実現した。

しかしSURFでもまた、特徴量の記述に高次元の実数を用いるためメモリ容量の増加と類似度計算の増加が問題となっていた。携帯端末向けの、より高速で省メモリな手法が求められ、特徴点検出のために、円周上の注目画素の輝度を決定木により学習するFAST、特徴量をバイナリコードで表現するBRIEF、BRISK、ORB、CARDが提案された。

その中でも、BRIEFを発展させた以下のような方法で、スケール不変性と回転不変性を得ているのがBRISKである。

ある画素の周りの円周上の画素を観測してエッジを検出するのがFASTであるが、BRISKでは入力画像を多段階に縮小したピラミッド画像のそれぞれに対して、FASTによる方法を適用して「特徴点らしさ」を求めている。この特徴点らしさがしきい値以上で、かつ、スケール変化方向において特徴点らしさが極大となる点を「特徴点」として抽出し、そのときのスケールを「特徴量」として求めている。

f:id:berobemin2:20151013163416p:plain

BRISKは、上図のようなパッチ(ある注目点に対して輝度比較を行う点を表わすもの)に配置された4つの同心円上において、等間隔にサンプリングされた全60箇所の画素値を用いている。

BRISKは、規則的な60個の決められた位置からランダムに2画素を選び、あらかじめスケールに応じて定められたしきい値をもとに「長距離ペアL」と「短距離ペアS」に分類し、それぞれのペアにおいて勾配に相当する量を計算する。パッチ内の大局的な勾配方向を捉えるために、長距離ペアLにおいて平均勾配ベクトルを求め、その角度を特徴点の「オリエンテーションα」とする。推定したオリエンテーションにしたがって輝度値を回転させて規格化する。規格化した輝度に対して、短距離ペアSの輝度比較を行い、バイナリコードの各ビットを算出する(1番目のペアの輝度比較結果0or1をi番目のバイナリに保存)。

これに対してBRIEFは、60個の画素からランダムに観測点を選んでいて、(ビット数×2)回分の画素値アクセスが必要となる。

なお、各サンプリング点における輝度は、中心からの距離に比例する分散をもつガウスフィルタによりスムージングをかけることでノイズ耐性を獲得している。

参考:http://www.vision.cs.chubu.ac.jp/CV-R/jpdf/Leutenegger_iccv2011.pdf

Cascaded FAST

以前の日記で、FASTではコーナー検出のために観測する画素を「決定木」の手法を用いて学習させていくと述べた。これを異なるサイズを組み合わせる方法に発展させたものが、Cascaded FASTと呼ばれる。

異なる周囲長を参照する3種類の決定木をカスケード状に並べて、学習により不要なコーナー点の検出を抑制し、コーナーらしい点のみを高速に検出する方法。オリエンテーションを算出してコーナーらしさを評価し、ピラミッド処理を用いてスケールを求める。このオリエンテーションとスケールを特徴量とし、回転不変性のある特徴点マッチングを行うことができる。

参照:http://mprg.jp/research/keypoint_j