lgammaの差を早く計算したい

ベイズを使った研究ではよく多項分布とディリクレ分布を使い、パラメータを積分消去したりするので、ポリヤ分布が頻出する。ポリヤ分布のcomponentはガンマ関数からなっており、普通は対数を取って計算することが多いので、ガンマ関数の対数の差(lgammaの差と以後書く)をよく計算することになる。collapsed LDAのようなケースでは対数尤度を計算するときだけこのlgammaの差が登場するが(sampling時には特に関係ないということ)、普通はここでボトルネックになることは少ない。しかし、例えばnested Chinese Restaurant Process(nCRP)のようなケースではパスのサンプリング式の一部で
$p(\mathbf{w}_d | \mathbf{c}_{d}, \mathbf{c}_{- d}, \mathbf{w}_{-d}, \mathbf{z}) = \prod_{l=1}^L \left(\frac{\Gamma{(\{N_{c_{d, l}}\}_{-d} + V \eta_l )}}{\prod_{v=1}^V \Gamma{(\{N_{c_{d, l}, v}\}_{-d} + \eta_l )}} \frac{\prod_{v=1}^V \Gamma{(N_{c_{d, l}, v} + \{N_{c_{d, l}, v}\}_{-d} + \eta_l)}} {\Gamma{(N_{c_{d, l}} + \{N_{c_{d, l}}\}_{-d} + V \eta_l)}} \right)$
というような式が登場し、サンプリングの度にlgammaの比を計算しなければならない(サンプリングのlogsumexpを使う関係で対数を取っている)。サンプリングの計算はMCMCで時間がかかるところなので、lgammaの差を高速にできるようになることは重要(なはず)。サンプリングの式にlgammaの差が登場するのはnCRPだけかというとそうではなくて単語レベル以外に例えば文に隠れ変数がいるようなモデルではサンプリング式にガンマ関数が登場するので、そこそこよくあるケースだと思う(ベイズ系修士学生なら在学中に3回くらいは遭遇するはず)。

そもそもnCRPのような例でなぜlgammaの差が登場してしまうかというと、ガンマ関数の性質
$\Gamma(n + 1) = n \Gamma(n)$
がうまく利用できないからである(キャンセルアウトでガンマ関数が消えてくれる、ということがない)。しかし、上の性質が全く利用できないかと言われるとそんなこともない。両辺に対数を取って移項すると
$\log \Gamma(n + 1) - \log \Gamma(n) = \log n$
という関係式を得ることができる。計算したいのはlgammaの差であるが、これを
$\log{\Gamma{(a + b)}} - \log{\Gamma{(a)}}$
とする。 $\log \Gamma(n + 1) - \log \Gamma(n) = \log n$ をうまく利用して、 $\log{\Gamma{(a + b)}} - \log{\Gamma{(a)}}$ を捻り出すことを考える。持ってる関係式は一個づつの差なのだからb回くらい*1繰り返せば出てきそうなので繰り返してみる。

$\log \Gamma(a + b) - \log \Gamma(a + b - 1) = \log{(a + b - 1)}$
$\log \Gamma(a + b - 1) - \log \Gamma(a + b - 2) = \log{(a + b - 2)}$
...
$\log \Gamma(a + 1) - \log \Gamma(a) = \log{(a)}$

bはそもそも整数でいいのかという疑問があるが、nCRPのようなケースだとbはXXXの回数となっているようなもので置けるので問題ない(ハイパーパラメータのほうは実数だが、それはaに入れてあげればよい)。上の式を全部足してみると、あれやこれやとキャンセルアウトが起こって
$\log{\Gamma{(a + b)}} - \log{\Gamma{(a)}} = \sum_{i=0}^{b-1} \log{(a + i)}$
という式を得ることができる。問題はこの対数の和がlgammaの差より早いのかどうかということであるが、手元で簡単に実験してみたところ対数の和にしたほうが3倍くらい早かった。実際の値でも一致していたので、上の計算もたぶん合ってるはず...。論文とか読んでると実際の実装どうしたとかあまり書いてなくて、こういう早くする話がどっかにあってもいいんじゃないかと思ったので書いてみた次第です。