tsujimotter

[latexpage]

＊

大きさが極端に小さい／大きい「重み」の値の和を求める際に、アンダーフロー／オーバーフローを防ぐための方法です。ベイズで周辺確率を求めるときなど計算機統計の分野でしばしば用いられます。

＊

応用の幅は広いと思いますが、今回はパーティクルフィルタという手法を例にとり、説明します。

＊

ここでパーティクルフィルタについての詳しい解説はしませんが、簡単に言うと、パーティクルフィルタは、重みのついたパーティクルと呼ばれる粒子を多数用意して、そのパーティクルの分布を使って任意の確率分布を近似する手法です。モンテカルロ法から出発しているので、モンテカルロフィルタとか逐次モンテカルロ法などと呼ばれることもあります。

パーティクルは最初は一様な分布で全面にばら撒かれます。その後、何らかの外的な基準により、分布を正しく表している（尤もらしい）と考えられる粒子の重みは大きくなり、逆に分布を正しく表していない（尤もらしくない）粒子はその重みが小さくなります。この試行を繰り返すことで、パーティクルの近似分布がより正確になっていくのです。この尤もらしさを表す指標は尤度と呼ばれています。

＊

さてこのパーティクルフィルタでは、重みの片寄りが大きくなり過ぎた場合に、その重みの大きさにしたがってパーティクルを確率的に再選択するリサンプリングという過程が存在します。

＊

具体的に、$i$番目のパーティクルの重みを$w_i$としましょう。この$i$番目のパーティクルが再選択される確率$P_i$は次式で表されます。

\begin{equation}

P_i = \frac{w_i}{\sum w_i}

\end{equation}

＊

さあ、ここからが問題です。

このパーティクルの重みは「非常に小さい」ものとします。

プログラミング言語には、倍精度実数型(double型)がありますが、その下限$10^{-308}$を下回ってしまえば計算上は$0$と変わりません。これをアンダーフローと言います。したがって、もしすべての重みがアンダーフローしてしまうと、$Q_i$の分母である重みの和は$0$となってしまい、$0$割りでエラーとなってしまいます。

全部が全部極端に小さくなかったとしても問題は発生します。倍精度実数の有効数字は15桁程度ですので、大きな重み和に小さな重みを足したときに、小さな重みが無視されてしまいます。

＊

そこで、重み$w_i$を対数として保持しておき、可能な限り対数のまま計算するという戦略がとられるわけです。

対数重みを$l_i=\log{w_i}$としたとき、$P_i$は次のように置き換えられます。

\begin{equation}

P_i = \frac{\exp{\log{w_i}}}{\sum \exp{\log{w_i}}}=\frac{\exp{l_i}}{\sum \exp{l_i}}

\end{equation}

最終的な結果も、対数で持っておくと便利なので、両辺対数をとります。

\begin{equation}

\log{P_i} = \log{\frac{\exp{l_i}}{\sum \exp{l_i}}}

\end{equation}

\begin{equation}

\log{P_i} =l_i-\log{\sum\exp{l_i}}

\end{equation}

対数をとった確率$P_i$が計算できました。ずいぶんシンプルな和になりましたね。

＊

これで、万事解決（？）と思いきや、そう簡単にはいきません。

たしかに、重みを対数で持つことにしたので、変数$l_i$の値が$0$である心配はありません。

しかしながら、式(4)の第二項目で$\exp{l_i}$の和を計算しており、和のアンダーフローの問題や有効数字の問題が残っています。これはなんとかしたいです。

＊

＊

実は、これはlogsumexpの問題と呼ばれており、計算機統計学の実にさまざまな分野で顔を出します。

式の形が$\log\sum\exp$と並んでいるために、logsumexpと呼ばれています。名前がついていることからわかるように、既にいろいろな研究者がこの問題に取り組んでおり、さまざなな解決方法が提案されています。

＊

＊

解決策１

解決策の１つとして、最小値を用いる方法があります。

ここで、$l_i$の最小値を$l_{min}$とします。

\begin{equation}

\l_{min}=\min {\left(l_1, l_2, \cdots, l_n\right)}

\end{equation}

これを使って、式(4)の第二項目を以下のように書き換えます。

\begin{eqnarray}

\log{\sum\exp{l_i}}&=&\log{\sum\exp{\left(l_i-l_{min}+l_{min}\right)}}\\

&=&\log{\sum\left(\exp{l_{min}}\exp{\left(l_i-l_{min}\right)}\right)}\\

&=&\log{\left[\left(\exp{l_{min}}\right)\left(\sum\exp{\left(l_i-l_{min}\right)}\right)\right]}\\

&=&l_{min}+\log{\sum\exp{\left(l_i-l_{min}\right)}}\\

\end{eqnarray}

式(4)に代入すると、

\begin{eqnarray}

\log{P_i}=\left(l_i-l_{min}\right)-\log{\sum\exp{\left(l_i-l_{min}\right)}}\\

\end{eqnarray}

式の形をよく見ると、これは式(4)の$l_i$を下記のように変数変換したものと同じです。

\begin{eqnarray}

l_i\rightarrow l_i-l_{min}

\end{eqnarray}

$l_i-l_{min}$の大きさが、$0$から$l_{max}-l_{min}$の間に収まるので、指数をとってもアンダーフローは起きません。

＊

＊

解決策２

もう一つのやり方は、二項間の和で関数を定義してしまう方法です。この方法が一般的かもしれません。

＊

${\rm logsumexp}(x, y)$という二項間関数を作って下記のように定義します。

\begin{eqnarray}

{\rm logsumexp}(x, y) = \log \left(\exp{(x)}+\exp{(y)}\right)

\end{eqnarray}

＊

指数の和はアンダーフローの問題がありますので、$x>y$を仮定し、次の変形を行います。

\begin{eqnarray}

{\rm logsumexp}(x, y) &=& \log \left(\exp{(x)}+\exp{(y)}\right)\\

&=& \log \left[\exp(x) \left(1+\exp{(y-x)}\right)\right]

\end{eqnarray}

これによって、指数部の計算がexp(y-x)に収まります。

＊

次のJavaコードで示す${\rm logsumexp}(x, y)$関数は、上記の要領で、logsumexpを計算する関数です。

public static double logsumexp (double x, double y, boolean flg)　{
    if (flg)　　 return y; 　　　　　　　　 // init mode
    if (x == y)　return x + 0.69314718055;　// log(2)
    double vmin = (x<=y) ? x : y;
    double vmax = (x>=y) ? x : y;
    if (vmax > vmin + 50) {
        return vmax;
    } else {
        return vmax + Math.log (Math.exp (vmin - vmax) + 1.0);
    }
}

（完全なソースコードはこちらStatistics.java）

＊

下式のような三項間の和を求める際も、

\begin{eqnarray}

\log (\exp(x) + \exp(y) + \exp(z))

\end{eqnarray}

二項間の和${\rm logsumexp}(x, y)$を用いて再帰的に計算可能です。

\begin{eqnarray}

{\rm logsumexp}\left({\rm logsumexp}(x, y) , z\right)&=&{\rm logsumexp}\left(\log \left(\exp{x}+\exp{y}\right) , z\right)\\

&=&\log \left(\exp{\log \left(\exp{x}+\exp{y}\right) }+\exp{z}\right)\\

&=&\log \left(\exp{x}+\exp{y}+\exp{z}\right)\\

\end{eqnarray}

＊

$n$ 個の指数和を For 文で求める場合も同様ですが、最初の和だけ注意する必要があります。

double[] x = {x1, x2, x3};
double sum;
for(int i=0; i<x.length; i++)
    sum = logsumexp(sum, x[i], (i==0));

＊

$i=0$ のときの和は通常、 $w$ が初期化されていないため不定になってしまいます。 logsumexp関数では、この問題を防ぐため、flg引数を用意し、初期状態のときのみ $y$ をそのまま返すようにしています。

＊

この方法は、下記のURLで紹介されていたものですが、

http://d.hatena.ne.jp/echizen_tm/20100628/1277735444

リンク先で述べられているように、もともとは、下記のPDF文献の後半に紹介されていたものです。

Sequential Tagging メモ 浅原正幸

リンクがなくなったときのために、コード自体も本記事に貼り付けさせて頂きました。

このような内容を紹介してくださっている皆様方に感謝です。

＊

＊

実はこの方法、計算機統計学の分野では浸透されていて常識になっているようですが、書籍や論文にはほとんど書かれていないのです。

実装のためのノウハウなので書籍や論文になりにくいのかもしれません。

統計屋さんの方々は研究室代々で伝承されていくのかもしれませんが、必ずしもそうでない場合もあるでしょう。多くの場合は、何とか試行錯誤して自力で見つけるか、あるいは不正確を承知で使っています。

実際私も周りに詳しい人がいないので、こんな方法があるという存在すらしりませんでした。

$0$割りのエラーが頻発して、どうしたものかと悩んでいた時に、たまたま持っていた計算機統計の本のあるページに「logsumexpという方法がある」と書いていたのです。しかも欄外に一行だけ。ネット上の情報をかき集め、なんとか実装しました。書かれている記事があまりに少ないのです。

＊

logsumexpの解決法にもまだまだ問題はあります。それは、$\log$や$\exp$の計算が頻発する点です。$\log$や$\exp$の計算速度は一般に遅いのです。

速度問題の解消のために「スケーリング」と呼ばれる方法も提案されているそうです。調べてみてください。私に説明できる力はありませんが、名前さえわかっていれば調べることは可能でしょう。

＊

車輪の再生産はもったいないと思いませんか。

＊

3 Responses to “logsumexp”