exp(x)の高速計算～理論編～

2010.09.02

exp(x)の高速計算～理論編～

ロジスティック回帰やCRFなどの対数線形モデルの学習でよく出てくるのが，expの計算です．これをSSE2を使って高速化するのが，今回のテーマです．まずは，背景の理論を説明します．

まず，指数関数 $e^x$ を２の指数関数 $2^a$ に変換することを考えます（なぜ２の指数関数かはいずれ分かります）．

$e^x = 2^a$

両辺の自然対数をとり， $a$ について解くと， $a = x/\log 2$ ．IEEE754など，２を基数とした指数部を採用している浮動小数点形式では，整数 $n$ に対して $2^n$ を容易に構築できるので，上式の実数解 $a$ の代わりに整数， $n = \left\lfloor x/\log 2 \right\rfloor$ を用い， $e^x$ の大まかな値を計算することを考えます．ただし， $\lfloor a \rfloor$ は $a$ を超えない最大の整数を表します．

さて， $a$ を $n$ で近似したときの誤差 $(a - n)$ の範囲は， $0 \leq a - n < 1$ ．誤差 $(a - n)$ に $\log 2$ を乗じたものを， $b$ と定義すると，

$b = (a - n) \log 2 = x - n \log2$

上式の両辺の指数をとると，

$e^x = 2^n e^b$

ここで， $n$ は $a$ を超えない最大の整数なので， $b$ の値域は， $0 \leq b < \log 2$ ．

これらのことから， $e^x$ は以下のステップで計算出来ることが分かります．

$n = \left\lfloor x/\log 2 \right\rfloor$ を計算する．
$b = x - n \log2$ を計算する．
$e^b$ を定義域 $(0 \leq b < \log 2)$ において，テーラー展開もしくは最良近似多項式などで近似計算する．
浮動小数点形式を利用して $2^n$ を構築し， $e^b$ との積を計算する．

ステップ1は簡単そうに見えますが，負の数を整数にキャストすると，0に近づいてしまうことに注意が必要です．つまり，doubleからintへのキャストは，正の値は切り下げ，負の値は切り上げになってしまいます．このことに注意すると，このステップは次のように書けます（a < 0で，ifによる分岐を作らないのが賢いやり方）．

a = LOG2E * x;
a -= (a < 0);
n = (int)a;

ステップ２はそのまま書けばOKです．

ステップ３では， $e^b$ を多項式に近似します．多項式近似というとテーラー展開が一般的ですが，「･･･点の周りについてテーラー展開すると」という表現の通り，ある点の近傍の近似は良いのですが，そこから離れた点では近似精度があんまりよくないという欠点があります．Cephesというライブラリでは，パデ近似（Padé approximant）を採用していますが，除算を１回だけ行う必要があります．最近のCPUは速くなったとはいえ，SSE2で倍精度のパック除算（DIVPD）をやるのに必要なレイテンシ・スループットは，70・70で，倍精度のパック乗算（MULPD）の7・6と比べると，除算が圧倒的に遅いことになります．そこで，今回は最良近似式を用います．なにやら難しそうに見えますが，Sollya というソフトウェアを使えば，簡単に最良多項式が求まります．以下の例では， $[0, \log 2]$ の範囲で，５次の最良多項式を求めています．

> remez(exp(x), 5, [0, log(2)]);
0.99999989311082729779536722205742989232069120354073 + x *
(1.00001092396453942157124178508842412412025643386873 + x *
(0.49981934577169208735732248650232562589934399402426 + x *
(0.16775408658617866431779970932853611481292418818223 + x *
(3.87412011356070379615759057344100690905653320886699e-2 + x *
1.185268231308989403584147407056378360798378534739e-2))))

多項式がホーナー形式（horner scheme）で出てくるので，至れり尽くせりです．これをそのままC言語に書き直せば，ステップ３の実装は完了です．

y = 1.185268231308989403584147407056378360798378534739e-2;
y *= b;
y += 3.87412011356070379615759057344100690905653320886699e-2;
y *= b;
y += 0.16775408658617866431779970932853611481292418818223;
y *= b;
y += 0.49981934577169208735732248650232562589934399402426;
y *= b;
y += 1.00001092396453942157124178508842412412025643386873;
y *= b;
y += 0.99999989311082729779536722205742989232069120354073;

最後のステップ４ですが，もちろんpow(2, n)なんてやってしまったら，せっかくの高速化が台無しです．浮動小数点形式IEEE754では，２を基数として指数部が構成されているので， $2^n$ はちょっとした工夫で一瞬で作れます．あんまり移植性がないコードですが， $2^n$ は次のように作ります．

typedef union {
    double d;
    unsigned short s[4];
} ieee754;

ieee754 u;
u.d = 0;
u.s[3] = (unsigned short)(((n + 1023) << 4) & 0x7FF0);

こうすると，u.dで $2^n$ の値にアクセスできます．なぜそうなるかは，IEEE754の仕様と睨めっこすればすぐに分かると思います．あとは，u.dとyの積を計算すれば， $e^x$ の計算が出来たことになります．

Não Aqui!