logsumexp(log sum exponential)とは

計算機科学

uchiumiさんが最近"log sum exponential"という単語を連呼していた。
「ろぐ さむ えくすぽおねんしゃる？なにそれおいしいの？」
状態だったのだが、最近LDA(Latent Dirichlet Allocation)を実装した際に、この"log sum exponential"問題に遭遇したのでメモしておく。


例えばある変数集合{Ai}があったときに

Pi = Ai / ΣAi

として最大1になるようにして確率化したいことがある。ここでPiやAiが非常に大きい、または小さい値を取る場合PiやAiをそのまま保持せずにlog(Pi)やlog(Ai)として持っておきたいことがある。
以下、lnPi = log(Pi)、lnAi = log(Ai)として扱う。

lnPi = lnAi - log(ΣAi)
     = lnAi - log(Σexp(lnAi)))

とすればPiとaiを常に対数の形でもっておけるので大きな値や小さなをとっても安心だ。
ここで2行目の第2項はどうやって計算すればいいのだろうか？ふつうにexp(lnAi)=Aiと展開してしまうとAiはそのまま持てないという扱いだったので困ってしまう。
実はこの第二項が冒頭で述べた"log sum exponential"というもの。logΣexp(x)という形をしているので"log sum exponential"だ。この部分の計算方法は浅原さんの資料

Sequential Tagging メモ 浅原正幸

が詳しい。というか実装を書いてくださっている。この実装ではlogsumexp()という名前の関数になっている。今回はこれについて説明を書いておく。
logsumexp()のポイントは以下のとおり。変数a,bがあってa < bという大小関係があったときに

exp(a) + exp(b) = exp(b){exp(a-b) + exp(0)}
                = exp(b){exp(a-b) + 1}

という関係があることを利用する。この関係は左辺の式から単純にexp(b)を共通項として括り出しているだけで

x^3 + x^5 = x^5 * {x^(-2) + x^0}

などと同じ。実際に簡単な数を入れるとわかりやすい。
話を戻すと、↑のexp(a) + exp(b)の式の対数をとると

    log{exp(a) + exp(b)}
  = log[exp(b){exp(a-b) + exp(0)}]
  = log{exp(b)} + log{exp(a-b) + 1}
  = b + log{exp(a-b) + 1}

となる。この最後の行の式にはexp(a)やexp(b)は出現せず、exp(a-b)のみが出現しているが、(a-b)はa、bより絶対値が小さいのでそのままexp(a)やexp(b)を計算するより現実的な値になる。
ここで

z = log{exp(a) + exp(b)}

とおくと

   log{exp(z) + exp(c)}
 = log{exp(log{exp(a) + exp(b)}) + exp(c)}
 = log{exp(a) + exp(b) + exp(c)}

となる。ここでa,b,cをlnAiと考えれば

   log{exp(a) + exp(b) + exp(c) + ...}
 = log{exp(lnA0) + exp(lnA1) + exp(lnA2) + ...}
 = logΣexp(lnAi)

となり無事に"log sum exponential"が計算できたことがわかる。
というわけで見たとおり普通に計算するより重い感じの処理になるので大きいデータだと使うのに抵抗がある。しかし、そもそも"log sum exponential"が必要になるような状況というのはデータが大きい時だ。うーん。困った。

※本記事の内容については、理解にあたって大部分uchiumiさんのお世話になった。thanks!