« | »

2010.09.02

exp(x)の高速計算 ~実験編~

最後に,e^xの計算がどのくらい高速化できたのか,実験してみました.アルゴリズムとして,以下のものを比較しました.

  • libcのexp関数
  • パデ近似を用いたもの(cephesの実装に近いです)
  • 5次~13次でテーラー展開を0を中心として行ったもの
  • 5次~13次の最良多項式近似を[-0.5, +0.5]の範囲で行ったもの
  • 5次~13次の最良多項式近似を[0, \log 2]の範囲で行ったもの
  • 5次~11次の最良多項式近似を[0, \log 2]の範囲で行ったものをSSE2で実装したもの

SSE2の実装で11次で打ち止めしているのは,11次と13次で演算精度の差が殆どないことが分かり,単に実装するのが面倒くさかったからだけです.0を中心としたテーラー展開や,最良多項式を[-0.5, +0.5]の範囲で構成する手法は,nの値を切り下げではなく,四捨五入(つまり0.5を足してから切り下げる)で求めます(Wapitisse_mathfun.hのexp計算ルーチンは四捨五入を使っています).

実験では,0を平均としたガウス分布で乱数を10,000,000個発生させ,その値のexpを全て計算するのに要する時間と,libcの演算結果を正しいものと仮定したときの誤差を測定しました.実験環境は,Intel(R) Xeon(R) CPU 5140 (2.33GHz) 上でDebian Linuxを動作させたアプリケーション・サーバーです.gccのバージョンは4.1.2で,コンパイルオプションは「-O3 -fomit-frame-pointer -msse2 -mfpmath=sse -ffast-math -lm」です.実験に用いたコードは,こちらからダウンロードもしくは閲覧できます.

横軸に計算時間,縦軸に誤差のRMS(二乗平均平方根)をプロットしたものを示します.Performance of fast exp computation

縦に一本線が入っているのが,libcの計算時間です.これを基準に誤差を測定しているので,エラーは0です.この線よりも左側に入っていないと,高速化の意味が無いことになります.各アルゴリズムは,左上が5次の多項式近似で,右下に向かって7次,9次,11次,13次と,精密かつ遅くなっていきます.だいたい,10^{-16}あたりで計算誤差の改善が止まっています.これは倍精度浮動小数点の有効桁数である16桁と大体一致しています.

計算精度に着目すると,テイラー展開 < 最良近似式 [-0.5, +0.5] < 最良近似式 [0, log 2] という順で,どの近似式も次数を上げていくと,10^{-16}付近まで誤差が改善します.ただ,アルゴリズムをCで実装してしまうと,libcと比較しても大した速度向上が得られません.5次まで近似精度を落としたとしても,2倍程度しか速くなりません.libcのexpの実装も相当に高速化されているので,一筋縄ではいかないようです.

これに対し,アルゴリズムをSSE2で実装したものは,4倍~8倍程度高速化されています.しかも,計算精度はSSE2化していないものとぴったり一致しているので,計算精度を犠牲にすることなく,高速化が達成できたことになります.e^xを大量に計算する状況では,11次の最良近似式を用いたSSE2ルーチンを使えば,演算精度を落とすことなく高速化できますし,近似式の精度が低くて良ければ,次数を下げてさらに高速化できます.次期バージョンのCRFsuiteでは,その他の改良も含めて,SSE2ルーチンを搭載する予定です.

Trackback URL

Comment & Trackback

No comments.

Comment feed

Comment





XHTML: You can use these tags:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>