二進対数の高速計算を用いた単語クラスタリングの高速化
先日実装した単語クラスタリングでは，相互情報量を計算する際の対数計算がボトルネックとなっている．``Premature optimization is the root of all evil'' とは良く言ったものだが，これ以上コードを更新するつもりもなかったので，高価な演算（対数計算）を単純な演算に置き換える定番の最適化を行ってみた．
今回の高速化の対象は，
float pmi (float p_ij, float p_i, float p_j) { return p_ij == 0 ? 0 : p_ij * log2 (p_ij / (p_i * p_j)); }
という関数．この手の計算を高速化するには，可能な引数に対する結果を事前に計算した lookup table を用意するのが常套手段*1．lookup table を用いた対数計算の高速化手法は，以下の論文で提案されている．
- A hardware-independent fast logarithm approximation with adjustable accuracy (IEEE International Symposium on Multimedia 2008)
肝は lookup table のサイズを小さくするため，浮動小数点数を指数部 (exponent) と仮数部 (mantissa) に分解し，仮数部から対数への写像を前計算して lookup table を構築するところ．
IEEE 754 では正の浮動小数点数 f = (1 + mantissa) * 2^(exponent - 127) の二進対数 log2 (f) は，log2 (f) = log2 (1 + mantissa) + exponent - 127 と表現されるので，以下のようにして仮数に対応するの二進対数を前計算する．
// union for manipulating exponent and mantissa in floating numbers union byte_4 { unsigned int u; float f; explicit byte_4 (unsigned int u_) : u (u_) {} byte_4 (float f_) : f (f_) {} }; const size_t PREC = 23; std::vector <float> LUT (1 << PREC, 0); byte_4 b (0x3f800000u); // 1 * 2^(127 - 127) = 1 for (size_t i = 0; i < (1 << PREC); ++i, b.u += (1 << (23 - PREC))) LUT[i] = log2 (b.f);
仮数は 23 bit なので，近似無しだと lookup table に (1 << 23) * sizeof (float) = 32M ほど必要だが，（PREC < 23 にして）末尾の bit を多少犠牲にすれば，lookup table のサイズは（CPU の）キャッシュメモリに載る程度に小さくできる．lookup table を使った対数計算は以下（簡易化のため，LUT と PREC は大域変数 or メンバ変数として省略）．
inline float log2_ (byte_4 b) { return LUT[(b.u & 0x7fffff) >> (23 - PREC)] + (b.u >> 23) - 127; }
この log2_ を使って，pmi を計算する関数を定義する．高コストな除算を（ついでに乗算も）展開した関数も定義してみた．
float pmi_1 (float p_ij, float p_i, float p_j) { return p_ij == 0 ? 0 : p_ij * log2_ (p_ij / (p_i * p_j)); } float pmi_2 (float p_ij, float p_i, float p_j) { return p_ij == 0 ? 0 : p_ij * (log2_ (p_ij) - log2_ (p_i * p_j)); } float pmi_3 (float p_ij, float p_i, float p_j) { return p_ij == 0 ? 0 : p_ij * (log2_ (p_ij) - log2_ (p_i) - log2_ (p_j)); }
この pmi_1, pmi_2, pmi_3 をそれぞれ先日実装した単語クラスタリングに組み込んで log2 を用いた実装と比べてみた．クラスタリングの入力となる n-gram は動的ダブル配列を使って Wikipedia のテキスト処理を高速化 - ny23の日記で使った Wikipedia 全文 (約1.53G) から列挙した．
|T| = 279,024,968; |V| = 1,594,517 baseline: MI=1.637430; 1804.9s (g++ 4.1.2), 1429.1s (icc 12.0.3) ------------------------------------------------------------------- K=64 pmi_1 | pmi_2 | pmi_3 prec MI time | MI time | MI time ------------------------|---------------------|-------------------- 12 1.551528 449.3s | 1.504427 413.0s | 1.443919 530.7s 13 1.570644 425.7s | 1.541163 410.5s | 1.489108 544.4s 14 1.587452 426.0s | 1.562740 425.7s | 1.519756 537.4s 15 1.602825 456.5s | 1.581760 462.1s | 1.546882 585.2s 16 1.615481 473.2s | 1.599146 483.7s | 1.569260 604.7s 17 1.623932 490.1s | 1.616639 504.9s | 1.592678 614.7s 18 1.629893 505.0s | 1.626043 511.2s | 1.605505 653.0s 19 1.633992 520.0s | 1.632256 555.0s | 1.620285 700.5s 20 1.635955 634.8s | 1.635385 734.1s | 1.629901 723.6s 21 1.636751 758.2s | 1.636549 915.4s | 1.634273 772.6s 22 1.637175 927.6s | 1.637081 1069.7s | 1.636427 831.8s 23 1.637333 1023.8s | 1.637329 1129.0s | 1.637043 832.6s ------------------------------------------------------------------- 64bit GNU/Linux with Intel X5560 2.80GHz CPU (8M Cache) * -ffast-math を最適化オプションに追加して再実験
というわけで，pmi_1 (PREC ~ 18) を用いれば log2 に基づく pmi を用いた場合に比べて，gcc をコンパイラとして用いた場合で約4倍弱，icc をコンパイラとして用いた場合でも約3倍弱ぐらい速くなった．PREC < 23とする場合は log2_ が近似になるので呼び出し回数は最小限にとどめたほうが良さそうだ．論文では，PREC = 14 辺りが推奨されていたが，このマシンではキャッシュメモリが大きいこともあり，PREC = 18 (末尾の 5 bit を捨てて，lookup table は 250K) ぐらいで十分高速となるようだ．
Taylor 展開や最良近似多項式，SSE などを併用すればさらに高速化できるかもしれないが，大変そうなのでこの辺で退散．泥臭いことに足を突っ込みたくない人は，Intel Math Kernel Library*2などを利用するのが良いかと．
なお，単純に速いクラスタリング手法が必要な人は，
- Web-Scale K-Means Clustering (WWW 2010; poster)
の実装 SofiaKMeans とかを用いると良いかも知れない．[追記] 一応試しておくか，と思って k-means をさらに速くする - ny23の日記 の k-means 実装と適当に比較してみたら，あまり処理速度は変わらなかった（あれ？）．というわけで，通常の k-means で繰り返しを早めに打ち切るのが一番良いかも．
