ハミング距離の計算はホントに速いのか？
これは@sakanazensen君が主催する『Computer Vision Advent Calendar 2013』の12/8の記事です。今年はあまり活発でないようなので、小ネタですが参戦しました。
はじめに
昨今のコンピュータビジョン・パターン認識分野で特徴ベクトルのバイナリベースの記述法が流行っています。その利点の一つとして、特徴ベクトル間の距離としてコンピュータにとって計算が容易な「ハミング距離」が使える、というものがあります。これはXOR演算と PopCount演算（いくつのビットが1かをカウントする演算）で構成されており、特に近年のCPUにはまず搭載されているベクトル計算命令セットの一つ「SSE4.2」の専用命令「POPCNT」が高速演算の根拠としてよく引き合いに出されます。二つともかなりプリミティブな命令ですから確かに高速に計算できそうな感じはします。しかしながら、例えばL2距離（ユークリッド距離）やL1距離（マンハッタン距離）も同様にベクトル計算は可能で、そのあたりとの比較はこれまで十分にされていないように思います。しかも小サイズな整数型であればそれなりの並列性は実現できるはずで、意外とL2・L1も凄く速いかもしれません。このことが原因で「ハミング距離の計算そのものが速いのか」あるいは「バイナリ化の次元削減によって計算が速いのか」の切り分けが曖昧になりモヤっと感が残ります。
というわけで本エントリでは、要素が整数型のベクトルについて、それらのL2距離・L1距離・ハミング距離のベクトル計算（SSEのintrinsic命令）における計算コストについて検証してみます。浮動小数点型ベースのベクトル間の内積演算については@sakanazensen君のブログが参考になります。こちらも是非。
ベンチマーク用のコードと簡単な解説
環境によっても結果は変わるはずなので、お手軽にテストできるようコードをGitHubにて public domain で公開しました。「Visual Studio 2010」向けですので「gcc」の場合は多少の変更が必要です*1。
コードは下記のmain.cppが全てになります。
コードの冒頭にある、変数Dがベクトルの次元数、変数Nが辞書ベクトルの数、ベクトルの要素はuchar型（１バイト）です。デフォルトでは D=128*2、N=400万としています。Nが大きすぎると落ちますので空気を読んで使ってください*3。ベクトルは乱数で生成し（少し時間がかかる）、一つのクエリベクトルに対して辞書中から最近傍ベクトルを全探索します。このときの距離計算をSSEあり／なしで比較します。使用したい距離は search()／search_sse() 内のコメントアウト部を切り替えることで変えれます。フェアに評価するために、各距離でベクトルのビット長を統一しています。つまりD=128の場合、L2・L1距離では128バイトベクトルとして、ハミング距離では128×8=1024ビットバイナリコードとして、それぞれ計算します。
ベンチマーク結果
ベンチマーク環境は以下の通り：
| CPU | Intel Core i5 M560 @ 2.67GHz x2 (Nehalem) |
|---|---|
| メモリ | 8GB |
| IDE | Visual Studio 2010 Professional |
| プラットフォーム | x64 |
| オプション | /O2 |
結果は以下のようになりました。【追記2013/12/09：ハミング距離のコードを見直しやや高速化しました。】
| 距離 | SSEなし | SSEあり |
|---|---|---|
| L2 (Euclidean) | 461ms | 127ms |
| L1 (Manhattan) | 576ms | 90ms |
| ハミング 32bits | |
|
| ハミング 64bits | |
予想に反して、L1距離が最も高速という結果になりました。次点で64bitベースのハミング距離でした。意外と有意な差があります。
考察
まず結果から言えることは、ハミング距離の計算そのものが特に高速というわけではないようです。この理由について、作っていて感じたことが２つあります。下記のL1距離とハミング距離の関数を見てください。
//L1距離計算 inline int dist_l1_sse(unsigned char* p,unsigned char* q) { __m128i t=_mm_setzero_si128(); for(unsigned d=0;d<D;d+=ALIGN) { t=_mm_add_epi32(t, _mm_sad_epu8( _mm_load_si128((__m128i*)(p+d)), _mm_load_si128((__m128i*)(q+d)) ) ); } return int(t.m128i_i64[0]+t.m128i_i64[1]); } //ハミング距離計算（64ビット版で x64 のみ対応） inline int dist_hamming64_sse(unsigned char* p,unsigned char* q) { long long result=0; unsigned long long* p64=reinterpret_cast<unsigned long long*>(p); unsigned long long* q64=reinterpret_cast<unsigned long long*>(q); for(unsigned d=0;d<D;d+=ALIGN) { // this is probably faster than using _mm_xor_si128(). result+=_mm_popcnt_u64((*p64++)^(*q64++)); result+=_mm_popcnt_u64((*p64++)^(*q64++)); } return int(result); }
第一にL1距離計算に用いる「_mm_sad_epu8()」が優秀すぎ、第二にハミング距離計算に用いる「_mm_popcnt_u64()」がややお粗末なこと、があります。前者のほうは差分をとって絶対値をとってそれらを一部加算してくれるという、かゆいところに手が届きまくりの命令になっています。本来なら差分計算で1bit符号が増えるぶんの対策が必要なのですが、そういうものも不要です。これはすごい。一方後者はなぜか64bit長までしか対応しておらず、SSEの128ビットレジスタから値をそのまま渡すということができません。数多くのモダンなアルゴリズムにとって PopCount は超重要命令であるにもかかわらずこの設計ぶりには何か理由があるのでしょうか。。。
以上のことから、バイナリ化の真の恩恵というのは「次元削減（コード長の削減）」にあるということになります。例えば今回のベンチマークは128×8=1024ビットにしましたが、実際には128ビット程度まで落としても精度が落ちないという報告もあります（もちろんデータにもよります）。仮に1/8サイズまでデータ量を削減できると、通信量の削減、ストレージの節約、メモリ上へのガン並べによるディスクアクセス回避、などなど多くの点でメリットが生じます。一般のベクトルでもPCAなどによって次元削減可能ですが、例えばPCA-SIFTでも128次元⇒32次元で1/4サイズどまりですし、ビット単位での長さ調整もできず比較的使い勝手がよくありません。また次元削減時の変換に要する計算コストでも、バイナリコードのほうに分があります。このあたりがバイナリ化の大きな利点でしょうか。
おわりに
昨日７時間かけて作った（L2がメンドクサかった）のですが、実際に手を動かして検証するのは大事だなと思った一日でした。L1がここまで速いというのは自分でも意外でした。
上で述べたように、絶対値誤差和命令のほうはかなりギリギリまで効率化されている一方、PopCount命令のほうはまだまだ改善の余地があります。今後のハードウェアの発展によってより効率的な命令が実現されれば、今回の実験をくつがえす結果となることも予想されます。楽しみに待つことにしましょう。
バグ報告やより効率的な実装アイデアなど遠慮なくコメントください。他の環境での実験報告も大歓迎です。
