$$
\newcommand{\b}[1]{\mathbf{#1}}
$$

ディリクレ分布の最尤推定

今回はディリ分布の最尤推定法による数値最適解を得る方法について説明したいと思う。
二項分布や指数分布、正規分布などの簡単な分布では解析解を得られが、ディレクレ分布のような高次元で複雑な確率分布になると最尤推定やMAP推定で解析解を求めることは基本的に出来ない。そこで、数値的な計算処理をして求める必要がある。

特にディレクレ分布はナイーブベイズや混合モデルなどの共役事前分布として活躍するので様々な分野で必要になる。
自分の場合は機械学習プロフェッショナルシリーズのトピックモデルの2章ユニグラムモデルを読んでいる時に必要になった。

そこでプロフェッショナルシリーズの書籍にも紹介されていた以下の論文の行間を埋めながら進めていく。

Estimating a Dirichlet distribution Thomas P. Minka 2000 (revised 2003, 2009, 2012)

また、この資料は**ディガンマ関数の漸近展開から基本的なことまで、なるべく省略せず一から全て説明している。**また、論文でこの文章に対応するのは2ページの途中までである。

定式化

ディレクレ分布は以下のような確率分布である。

$$
\begin{eqnarray}
P(\b{p} | \b{\alpha}) = \frac{\Gamma(\sum_{k=1}^n \alpha_k)}{\prod_{k=1}^{n} \Gamma(\alpha_k)} \prod_{k=1}^{n} p_k^{\alpha_k - 1} \
\sum_{k=1}^{n} p_k = 1, ~~~~ p_k \ge 0
\end{eqnarray}
$$

ここで $ \b{\alpha} $ はn次元ベクトルである。

最尤推定法で最適解を求める

さて $ \b{\alpha} $ を固定してディレクレ分布から以下のように N 個の確率ベクトル $ DataSet = \{ \b{p^{(1)}}, \b{p^{(2)}}, ..., \b{p^{(N)}} \} $ が観測出来た時にその確率を最大化させるような $ \b{\alpha} $ を求める。尤度を $ L(\b{\alpha}) $ とすると以下のように表せる。

$$
\begin{eqnarray}
L(\b{\alpha}) & = & log(\prod_{i=1}^{N} P(\b{p^{(i)}} | \b{\alpha})) \
& = & \sum_{k=1}^{N} \{~ (\sum_{i=1}^{N} (\alpha_k - 1) ~ log(p_k^{(i)})) + N ~ log( \Gamma(\sum_{k'=1}^{n} \alpha_{k'})) - N ~ log(\Gamma(\alpha_k)) ~\} \
& = & \sum_{k=1}^{N} \{ ~ N ~ log( \Gamma(\sum_{k'=1}^{n} \alpha_{k'})) + ~ N (\alpha_k - 1) ~ log(\tilde{p_k}) - N ~ log(\Gamma(\alpha_k)) ~\}
\end{eqnarray}
$$

ここで $ log(\tilde{p_k}) $ については以下のように定義した。

$$
\begin{eqnarray}
log(\tilde{p_k}) = \frac{1}{N} \sum_{i=1}^{N} log(p_k^{(i)})
\end{eqnarray}
$$

さて $ L(\alpha) $ の極値を得るために、まず関数を微分してみる。

$$
\begin{eqnarray}
\frac{\partial L(\b{\alpha})}{\partial \alpha_i} & = & N \sum_{k=1}^{n} \{ \frac{\partial \alpha_k}{\partial \alpha_i} log(\tilde{p_k}) - \frac{\partial log(\Gamma(\alpha_k))}{\partial \alpha_i} \} + N \frac{\partial}{\partial \alpha_i} \Gamma(\sum_{k=1}^{n} \alpha_k) \
& = & \sum_{k=1}^{n} \{ \delta_{k,i} log(\tilde{p_k}) - \delta_{k, i} \psi(\alpha_k) \} + N ~ \psi(\sum_{k=1}^{n} \alpha_k) \frac{\partial \sum_{k=1}^{n} \alpha_k}{\partial \alpha_i} \
& = & N \{ \psi(\sum_{k=1}^{n} \alpha_k) - \psi(\alpha_i) + log(\tilde{p_i}) \}
\end{eqnarray}
$$

ここで $ \phi(x) $ をLogスケールのガンマ関数の微分、ディガンマ関数として新しく定義した。

$$
\begin{eqnarray}
\phi(x) = \frac{d}{dx} log(\Gamma(x)) = \frac{\Gamma'(x)}{\Gamma(x)}
\end{eqnarray}
$$

$ L(\b{\alpha}) $ は凸性を持っていることを示せるので極値は一つ。
すなわち

$$
\begin{eqnarray}
\frac{\partial L(\b{\alpha})}{\partial \b{\alpha}} = \b{0}
\end{eqnarray}
$$

を満たす点が最大値を取る $\alpha$ になる。
解析解がないのでこの値を数値解析で求める。不動点反復法とニュートン法、最急勾配法の3つの方法が利用できる。論文では不動点反復法とニュートン法についての説明をしているのでそれに従う。

不動点反復法による解法

まず一般的な不動点反復法の1変数の場合について説明する。
この方法は以下の様な不等式を満たす x に対して

$$
x = f(x)
$$

次のような漸化式を繰り返し用いることで答えを求める方法である。

$$
x^{(t)} = f(x^{(t-1)})
$$

もちろんどのような関数 $ f(x) $ に対しても上記の漸化式を用いて答えを求めることができるわけではない。リプシッツ連続性という性質を満たす場合のみ可能であるがその証明はここでは行わない。
さてこの関係性を多次元変数の関数である (11) の方程式に対しても成り立つと仮定すると以下の式を用いて更新すれば目的の値を求められそうである。

$$
\begin{eqnarray}
\psi(\alpha_i^{(t)}) & = & \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) + log(\tilde{p_i}) \
\alpha_i^{(t)} & = & \psi^{-1}(\psi(\sum_{k=1}^n \alpha_k^{(t-1)}) + log(\tilde{p_i})) \
\end{eqnarray}
$$

これを示すために以下のセクションの Appendix B で与えた不等式を利用して尤度の下限を求めてその下限の最大の求めることで、上記の更新式と同じ式を得られることを示す。 $ \tilde{x} = \sum_{k=1}^n \alpha_k^{(t-1)} $ として固定した時、(5)の一番目の項に対し Appendix B の不等式を適用することで

$$
\begin{eqnarray}
\frac{L(\b{\alpha})}{N} & \ge & log \{ \Gamma(\sum_{k=1}^n \alpha_k^{(t-1)}) exp((\sum_{k=1}^n \alpha_k - \alpha_k^{(t-1)})\psi(\sum_{k=1}^n \alpha_k^{(t-1)}))) \} \nonumber \
& & - \sum_{k=1}^{n} log(\Gamma(\alpha_k)) + \sum_{k=1}^n (\alpha_k-1)log(\tilde{p_k}) \nonumber \
& = & (\sum_{k=1}^n \alpha_k) \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) - \sum_{k=1}^{n} log(\Gamma(\alpha_k)) + \sum_{k=1}^n (\alpha_k-1)log(\tilde{p_k}) + const(\alpha_k^{(t-1)})
\end{eqnarray}
$$

ここで $ \alpha_k^{(t-1)} $ は定数とみなしている。よって、これらの項のみの部分を $ const(\alpha_k^{(t-1)}) $ でまとめた。
上記の不等式の右辺を $ G(\b{\alpha}) $ として $ \alpha_1, \alpha_2, ... \alpha_k $ について最大化させると不動点反復法の更新式と同じ式を得る。

$$
\begin{eqnarray}
\frac{\partial}{\partial \alpha_i} G(\b{\alpha}) & = & \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) - \sum_{k=1}^n \frac{\partial \alpha_k}{\partial \alpha_i} \psi(\alpha_k) + \sum_{k=1}^n \frac{\partial \alpha_k}{\partial \alpha_i} log(\tilde{p_k}) \
& = & \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) - \psi(\alpha_i) + log(\tilde{p_i})
\end{eqnarray}
$$

微分 = 0 で停留点が得られるので

$$
\begin{eqnarray}
\frac{\partial}{\partial \alpha_i} G(\b{\alpha}) & = & 0 \
\psi(\sum_{k=1}^n \alpha_k^{(t-1)}) - \psi(\alpha_i) + log(\tilde{p_i}) & = & 0 \
\psi(\alpha_i) & = & \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) + log(\tilde{p_i})
\end{eqnarray}
$$

このように下限値の最大化でも不動点反復法を仮定して導出した更新式を同じ式が得られた。（下限値の最大化を繰り返し行うという意味でEMアルゴリズムに似ている）

ここで下限の $ G(\b{\alpha}) $ を最大値をとるパラメーターが $ \psi(\alpha_i) = \psi(\alpha_k^{(t)}) $ となっている。

ニュートン法による解法

ニュートン法による解法も紹介する。
多次元のニュートン法は以下を満たす $ \b{x} $ を得る方法である。

$$
\begin{eqnarray}
\b{f}(\b{x}) = \b{0}
\end{eqnarray}
$$

$ \b{f}(\b{x}) $ を適当な点 $ \b{x^{\ast}} $ に対してテイラー展開すると

$$
\begin{eqnarray}
\b{f}(\b{x}) = \b{f}(\b{x^{\ast}}) + \frac{\partial \b{f}(\b{x^{\ast}})}{\partial \b{x}}(\b{x} - \b{x^{\ast}}) + O((\b{x} - \b{x^{\ast}})^2)
\end{eqnarray}
$$

二次以降の微小量を無視して$\b{f}(\b{x}) = \b{0}$ を代入し、 $\b{x}$ について解くと以下の式を得る。

$$
\begin{eqnarray}
\b{x} = \b{x^{\ast}} - \left(\frac{\partial \b{f}(\b{x^{\ast}})}{\partial \b{x}}\right)^{-1} \b{f}(\b{x^{\ast}})
\end{eqnarray}
$$

これがニュートン法の更新式である。

さて今求めたい値は $ \frac{\partial L(\b{\alpha})}{\partial \b{\alpha}} = \b{0} $ なので、$ L(\b{\alpha}) $ をもう一度微分する。

$$
\begin{eqnarray}
\frac{1}{N} \frac{\partial^2 L(\b{\alpha})}{\partial \alpha_j \partial \alpha_i} & = & \psi'(\sum_{k=1}^n \alpha_k) (\sum_{k=1}^n \frac{\partial \alpha_k}{\partial \alpha_j}) + \psi'(\alpha_i) \frac{\partial \alpha_i}{\partial \alpha_j} \
& = & \psi'(\sum_{k=1}^n \alpha_k) - \psi'(\alpha_i) \delta_{i,j}
\end{eqnarray}
$$

$\{ \b{H(\b{\alpha})} \}_{i,j} = \frac{\partial L(\b{\alpha})}{\partial \alpha_j \partial \alpha_i}$ として新しくヘッセ行列を定義する。また、 $ z = N ~ \psi'(\sum_{k=1}^n \alpha_k) $ と $ \{ \b{Q} \}_{i,j} = - N ~ \psi'(\alpha_i) \delta_{i,j} $, $ \{\b{g}(\b{\alpha})\}_{i} = \frac{\partial L(\alpha)}{\partial \alpha_i} $ と新しい行列・関数・ベクトル $ z, \b{Q}, \b{g} $ を定める。ヘッセ行列 $ \b{H} $ は以下のように表せる。

$$
\begin{eqnarray}
\b{H} = \b{Q} + \b{1 1^t}z
\end{eqnarray}
$$

$ L(\b{\alpha}) $ に対してニュートン法を適用すると以下のように表せる。

$$
\begin{eqnarray}
\b{\alpha^{(t)}} = \b{\alpha^{(t-1)}} - \b{H}^{-1}(\b{\alpha^{(t-1)}}) \b{g}(\b{\alpha^{(t-1)}})
\end{eqnarray}
$$

更新式では $ \b{H}, \b{g} $ がそれぞれ $ \b{x^{(t-1)}} $ に依存していることを明示的に示すために括弧を入れた。ニュートン法では $ \b{H} $ の逆行列を用いるので、これを解析的に表したい。そのために Appendix C で説明している Sherman–Morrison の公式を用いて式の変形を行うと

$$
\begin{eqnarray}
\b{H}^{-1}\b{g} & = & (\b{Q} + \b{1 1^t} z)^{-1}\b{g} \
& = & \b{Q}^{-1}\b{g} - \frac{\b{Q^{-1} 1 1^t Q^{-1} g}}{\frac{1}{z}+\b{1^t Q^{-1} 1}} \
& = & \b{Q}^{-1}\b{g} - \b{Q^{-1} 1} (\frac{\b{1^t Q^{-1} g}}{\frac{1}{z}+\b{1^t Q^{-1} 1}}) \
& = & \b{Q}^{-1} (\b{g} - \b{1}b)
\end{eqnarray}
$$

ここで $ b $ は以下の定数である。

$$
\begin{eqnarray}
b = \frac{\b{1^t Q^{-1} g}}{\frac{1}{z}+\b{1^t Q^{-1} 1}}
\end{eqnarray}
$$

$ \b{H}^{-1} \b{g} $　の各項を具体的に求めると　$ \{\b{Q}\}_{i,j} = - N ~ \psi'(\alpha_i) \delta_{i,j} $ を直接代入することで以下のようになる。

$$
\begin{eqnarray}
\{\b{H}^{-1} \b{g}\}_{k} & = & \{ \b{Q}^{-1} (\b{g} - \b{1}b) \}_{k} \
& = & \sum_{j} \{\b{Q^{-1}}\}_{k,j} (\b{g} - \b{1} b)_{j} \
& = & \sum_{j} \frac{\delta_{k,j}}{-N ~ \psi'(\alpha_k)} (g_j - b) \
& = & \frac{g_k - b}{-N ~ \psi'(\alpha_k)} \
& = & \frac{g_k - b}{\{\b{Q}\}_{k,k}}
\end{eqnarray}
$$

ニュートン法まとめ

これらの結果をまとめるとニュートン法の更新式は以下のように簡潔に表せる。

$$
\begin{eqnarray}
\b{\alpha}^{(t)} & = & \b{\alpha}^{(t-1)} - \b{H^{-1} g} \
\{\b{H^{-1} g}\}_{k} & = & (g_k - b)\{\b{Q^{-1}}\}_{k,k} \
b & = & \frac{\b{1^t Q^{-1} g}}{\frac{1}{z} + \b{1^t Q^{-1} 1}} \
\{\b{g}\}_{k} & = & N ~ \psi(\sum_{k=1}^n \alpha_k^{(t-1)}) - N ~ \psi(\alpha_k^{(t-1)}) + N ~ log(\tilde{p_k}) \
z & = & N ~ \psi'(\sum_{k=1}^n \alpha_k^{(t-1)}) \
\{\b{Q}^{-1}\}_{i,j} & = & - \frac{\delta_{i,j}}{N ~ \psi'(\alpha_i^{(t-1)})} \
\end{eqnarray}
$$

Appendix List

Appendix A: ディガンマ関数の計算と性質

ディガンマ関数の漸近展開

最初はディガンマ関数の具体的な級数表示を導く。

ディガンマ関数の定義は $ \frac{\partial}{\partial x} log \Gamma(x) $ である。$ log $ が入っているので、まずガンマ関数を掛け算の乗積表示として表せると便利であろう。これはワイエルシュトラスの無限乗積表示を用いることで解決する。

ワイエルシュトラスの無限乗積表示は以下の関係式のことを指す

$$
\begin{eqnarray}
\Gamma(x) & = & \lim_{n \to \infty} \frac{n!n^x}{(n+x)(n+x-1)...(x+1)x}
\end{eqnarray}
$$

この式の導出方法についても、数学的な厳密性は犠牲にしてできるだけ直感的に説明する。まず階乗を表すガンマ関数が与えられているとしよう。**このガンマ関数の具体的な形を知らなかったとしてガンマ関数を自然数から実数に拡張したい。**そのようなケースとして考えよう。

最初の出発点は以下の近似式である。
$ k << n $ の時

$$
\begin{eqnarray}
\frac{\Gamma(n+k+1)}{\Gamma(n+1)} & = & (n+k)(n+k-1)...(n+1) \
& \cong & n^k
\end{eqnarray}
$$

である。ここで $ k $ に対して自然数から実数への拡張を試みる。上記の関係を実数 $ x $ についても成り立つと拡張するのは自然であろう。すなわち $ x << n $ として

$$
\begin{eqnarray}
\frac{\Gamma(n+x+1)}{\Gamma(n+1)} & \cong & n^x
\end{eqnarray}
$$

十分大きな n について 上記の式が成り立つと仮定し、それぞれ両辺を展開すると

$$
\begin{eqnarray}
\lim_{n \to \infty} \frac{(n+x)(n+x-1)...(x+1)x\Gamma(x)}{n!} & = & \lim_{n \to \infty} n^x \
\Gamma(x) & = & \lim_{n \to \infty} \frac{n!n^x}{(n+x)(n+x-1)...(x+1)x}
\end{eqnarray}
$$

とワイエルシュトラスの無限乗積表示を導くことができる。

さてこれを用いてディガンマ関数の具体的な形を求める。

$$
\begin{eqnarray}
\frac{\partial}{\partial x} log(\Gamma(x)) & = & \lim_{n \to \infty} \frac{\partial}{\partial x} log(\frac{n!n^x}{\prod_{i=0}^n (i+x)}) \
& = & \lim_{n \to \infty} \frac{\partial}{\partial x} \{
\sum_{i=1}^n log(i) - log(i+x) \} - log(x) + x ~ log(n) \
& = & \lim_{n \to \infty} \{- \sum_{i=1}^n \frac{1}{i+x} - \frac{1}{x} + log(n) \} \
& = & - \frac{1}{x} + \lim_{n \to \infty} \{ (\sum_{i=1}^n \frac{1}{i} - \frac{1}{i+x}) - (\sum_{i=1}^n -\frac{1}{i} + log(n)) \} \
& = & -\gamma - \frac{1}{x} + (\sum_{i=1}^{\infty} \frac{1}{i} - \frac{1}{i+x}) \
\end{eqnarray}
$$

ただし $ \gamma $ はオイラー定数で以下の極限である。

$$
\begin{eqnarray}
\gamma = \lim_{n \to \infty} (log(n) - \sum_{i=1}^n \frac{1}{i})
\end{eqnarray}
$$

ここで以下の関係式が成り立つことに注目すると

$$
\begin{eqnarray}
\phi(x + 1) & = & \frac{\partial}{\partial x} log \Gamma(x + 1) \
& = & \frac{\partial}{\partial x} (log(x) + \Gamma(x)) \
& = & \frac{1}{x} + \phi(x)
\end{eqnarray}
$$

先ほどの式は以下のように表せる。

$$
\begin{eqnarray}
\psi(x+1) = -\gamma + (\sum_{i=1}^{\infty} \frac{1}{i} - \frac{1}{i+x})
\end{eqnarray}
$$

このように具体的なディガンマ関数の級数表示が得られた。
また二階微分は以下のようになる。

$$
\begin{eqnarray}
\psi'(x+1) = \sum_{i=1}^{\infty} \frac{1}{(i+x)^2} \
\end{eqnarray}
$$

この二階微分は逆2乗和なので収束は極めて遅い。
そのためより収束が早い関数の級数として表したい。そのためにオイラーマクローリンの和公式を使ってベルヌーイ数の級数として表す。

オイラーマクローリンの和公式

オイラーの和公式を説明する。微分演算子を $ D = \frac{d}{dx} $ とする。
この時テイラー展開を用いると

$$
\begin{eqnarray}
f(x+n) & = & \sum_{k=0}^{\infty} \frac{D^k f(x)}{k !} n^k \
& = & (\sum_{k=0}^{\infty} \frac{(n D)^k }{n!}) f(x) \
& = & e^{nD} f(x)
\end{eqnarray}
$$

これより

$$
\begin{eqnarray}
\sum_{k=0}^{N-1} f(x+n) & = & \sum_{n=0}^{N-1} e^{nD} f(x) \
& = & \frac{e^{N ~ D} - 1}{e^D - 1} f(x) \
& = & \frac{1}{e^D - 1} (f(x + N) - f(x)) \
& = & (D^{-1} + \sum_{n=1}^{\infty} \frac{B_n}{n!} D^{(n-1)}) (f(x + N) - f(x))
\end{eqnarray}
$$

ここで以下のように母関数のテイラー展開を用いてベルヌーイ関数 $ B_n $ を定義した。

$$
\begin{eqnarray}
\frac{x}{e^x - 1} = \sum_{n=0}^{\infty} \frac{B_n}{n!} x^n
\end{eqnarray}
$$

両辺の $ \lim{x \to 0} $ の極限を考えれば $ B_0 = 1 $ であることが分かる。また $ D^{-1} $ の逆演算の積分であることは次のようにして理解できる。

$$
\begin{eqnarray}
g(x) & = & D^{-1} f(x) \
D g(x) & = & D \circ D^{-1} f(x) \
\frac{d}{dx} g(x) & = & f(x) \
g(x) & = & \int_{\alpha}^{x} f(t) dt + C
\end{eqnarray}
$$

これらより元の方程式に $ x = 0 $ を代入して N に対して極限を取ると

$$
\begin{eqnarray}
\sum_{k=0}^{\infty} f(n) - \int_{0}^{\infty} f(n) dn & = & \lim_{N \to \infty} C_N + \sum_{n=1}^{\infty} \frac{B_n}{n!} (f^{(n-1)}(N) - f^{(n-1)}(0))
\end{eqnarray}
$$

母関数同士の引き算を計算することで3以上の奇数で全て $ B_n = 0 $ になることを証明する。 $ g(x) = \frac{x}{e^x - 1} $ として

$$
\begin{eqnarray} \
g(x) - g(-x) & = & \frac{x}{e^x - 1} - \frac{-x}{e^{-x} - 1} \
& = & \frac{x(1-e^x)}{e^x-1} = -x
\end{eqnarray}
$$

これより

$$
\sum_{k=0}^{\infty} \frac{2 B_{2k+1}}{(2k+1)!} x^{2k} = -x
$$

で両辺を比べることで $ B_1 = - \frac{1}{2} $ かつ $ B_n = 0, n \ge 3 $ であることが分かる。これを用いてオイラーマクローリンの公式を書きなおせば以下のようになる。

$$
\begin{eqnarray}
\sum_{k=0}^{\infty} f(n) - \int_{0}^{\infty} f(n) dn & = & \lim_{N \to \infty} C_N + \frac{1}{2} (f(N) - f(0)) + \nonumber\
& & ~~~~~~ \sum_{n=1}^{\infty} \frac{B_{2n}}{(2n)!} (f^{(2n-1)}(N) - f^{(2n-1)}(0))
\end{eqnarray}
$$

ここで $ C $ は適当な定数である。（右辺の極限が収束する場合の剰余項）

ディガンマ関数の漸近表示

オイラーマクローリンの公式に対して以下のように $ f(n) = \frac{1}{(n+x)^2} $ を定義すればディガンマ関数の微分は以下のように表せる。

$$
\begin{eqnarray}
\psi'(x) - \int_{0}^{\infty} \frac{1}{(n+x)^2} dn & = & C + \frac{1}{2x^2} + \sum_{n=1}^{\infty} \frac{B_n}{x^{2n+1}} \
\psi'(x) & = & C + \frac{1}{x} + \frac{1}{2x^2} + \sum_{n=1}^{\infty} \frac{B_{2n}}{x^{2n+1}}
\end{eqnarray}
$$

両辺 $ \lim x \to \infty $ の極限を考えれば $ C = 0 $ であることが分かる。

$$
\begin{eqnarray}
\psi'(x) & = & \frac{1}{x} + \frac{1}{2x^2} + \sum_{n=1}^{\infty} \frac{B_{2n}}{x^{2n+1}}
\end{eqnarray}
$$

この両辺を積分すれば目的のディガンマ関数の級数表示が得られる。

$$
\begin{eqnarray}
\psi(x) & = & C + log(x) - \frac{1}{x} - \sum_{n=1}^{\infty} \frac{B_{2n}}{2n ~ x^{2n}}
\end{eqnarray}
$$

ここで $ \psi(x+1) = \psi(x) + \frac{1}{x} $ を繰り返し適用することで自然数 M に対して以下のように閉じた形で表すことが出来る。

$$
\begin{eqnarray}
\psi(M) & = & \sum_{k=1}^{M-1} \frac{1}{k} + \psi(1) \
& = &\sum_{k=1}^{M-1} \frac{1}{k} - \gamma
\end{eqnarray}
$$

したがってこのような自然数 M を先ほどの式に代入してその極限を取ると

$$
\begin{eqnarray}
\lim_{M \to \infty} \psi(M) - log(M) & = & C + \lim_{M \to \infty} - \frac{1}{M} - \sum_{n=1}^{\infty} \frac{B_{2n}}{2n ~ M^{2n}} \
C & = & \lim_{M \to \infty} -\gamma + \{ \sum_{k=1}^{M-1} \frac{1}{k} - log(M) \} \
& = & - \gamma + \gamma = 0 \
\end{eqnarray}
$$

よって定数は 0 になることが分かる。
これらをまとめるとディガンマ関数とディガンマ関数の微分 $ \psi(x), \psi'(x)$ はそれぞれ以下の級数を数値計算すればよい。またベルヌーイ数の求め方については Appendix D にて説明している。

$$
\begin{eqnarray}
\psi'(x) & = & \frac{1}{x} + \frac{1}{2x^2} + \sum_{n=1}^{\infty} \frac{B_{2n}}{x^{2n+1}} \
\psi(x) & = & log(x) - \frac{1}{x} - \sum_{n=1}^{\infty} \frac{B_{2n}}{2n ~ x^{2n}}
\end{eqnarray}
$$

Appendix B: ディガンマ関数の不等式

証明するべき不等式は

$$
\begin{eqnarray}
\Gamma(x) \ge \Gamma(\tilde{x})~exp((x-\tilde{x}) \psi(x))
\end{eqnarray}
$$

である。
この不等式を証明する前にまずディガンマ関数の単調性を利用する。この事実は　$ \psi'(x+1) = \sum_{i=1}^{\infty} \frac{1}{(i+x)^2} > 0 $ から明らかである。この単調増加性を平均値の定理と組み合わせれば目的の不等式が得られる。

任意の $ y > x $ に対して平均値の定理より

$$
\begin{eqnarray}
\frac{log(\Gamma(y)) - log(\Gamma(x))}{y-x} = \left.\frac{d~ log(\Gamma(x))}{dx}\right|_{x=z} = \psi(z) ~~~~ (x < z < y)
\end{eqnarray}
$$

ここでディガンマ関数の単調性を利用すると

$$
\begin{eqnarray}
\frac{log(\Gamma(y)) - log(\Gamma(x))}{y-x} > \psi(x)
\end{eqnarray}
$$

同様にして $ y < x $ に対しても

$$
\begin{eqnarray}
\frac{log(\Gamma(x)) - log(\Gamma(y))}{x-y} = \left.\frac{d~ log(\Gamma(x))}{dx}\right|_{x=z} = \psi(z) ~~~~ (y < z < x)
\end{eqnarray}
$$

より単調性から以下が成り立つ

$$
\begin{eqnarray}
\frac{log(\Gamma(x)) - log(\Gamma(y))}{x-y} < \psi(x)
\end{eqnarray}
$$

これらの結果を合わせれば以下の不等式が任意の $ \tilde{x} $ に対して成り立つことは明らかである。

$$
\begin{eqnarray}
log(\Gamma(x)) - log(\Gamma(\tilde{x})) \ge \psi(x)(x-\tilde{x})
\end{eqnarray}
$$

両辺の指数をとって式を変形させると目的の不等式を得る。

$$
\begin{eqnarray}
\Gamma(x) \ge \Gamma(\tilde{x})~exp((x-\tilde{x}) \psi(x))
\end{eqnarray}
$$

Appendix C: Sherman–Morrison の公式導出

Sherman Morrison の公式とは以下のような関係式のことである。

$$
\begin{eqnarray}
(\b{A} + \b{u} \b{v^t})^{-1} & = & \b{A}^{-1} - \frac{\b{A}^{-1}\b{u}\b{v^t}\b{A^{-1}}}{1+\b{v^t}\b{A^{-1}}\b{u}}
\end{eqnarray}
$$

左辺の形をした式で $ \b{A^{-1}} $ の逆行列の値が既に分かっている場合、左の式ではなく右の式で計算をすると逆行列を求めることなく目的の値を得ることができる。なぜならば、逆行列はガウスジョルダンの消去法やLU分解な基本的に $ O(n^3) $ で重い処理なのでなるべく計算したくないからである。また扱うパラメーターが多い場合そもそも求めることが出来ない。

さて、上記の式を証明するために天下り的に $ \b{A} + \b{u}\b{v^t} $ を掛けて恒等行列になることを確認してもよいのだが、そのような方法では必要な時に自分で導出できないので困る。

そこで以下ではノイマン級数を用いた方法と代数的な方程式に帰着させて求める方法の二つを紹介する。自分のオススメは（厳密性は多少欠けるものの）ノイマン級数を使った方法である。

この方法さえ覚えておけば woodbury の公式など逆関数で表せれた式の変形が自由にできるようになる。その上、未知の問題に出会った時も応用がきくからである。

ノイマン級数を用いた方法

ノイマン級数とは以下のことで

$$
\begin{eqnarray}
(\b{A} - \b{I})^{-1} & = & \b{I} + \b{A} + \b{A^{2}} + \b{A^{3}} + ....
\end{eqnarray}
$$

$ |x| < 1 $ において成り立つ無限級数を行列で表したものである。

$$
\begin{eqnarray}
\frac{1}{1-x} & = & 1 + x^2 + x^3 + ....
\end{eqnarray}
$$

ノイマン級数も $ |\b{A}| < 1 $ が成り立たなければならず、成り立たない場合は右辺の無限級数が発散するためノイマン級数の等式を使用してはならない。 なのでこの等式を途中で使う場合、この制約を満たしている行列にしか証明したことにならないことは注意するべきである。

さて証明自体は簡単で以下のように式変形を行えばよい。

$$
\begin{eqnarray}
(\b{A} + \b{u} \b{v^t} )^{-1} & = & \{ \b{A} (\b{I} + \b{A^{-1} \b{u} \b{v^t}}) \}^{-1} \
& = & (\b{I} + \b{A^{-1}}\b{u}\b{v^t}) \b{A}^{-1} \
& = & (\sum_{k=0}^{\infty} (-\b{A^{-1}}\b{u}\b{v^t})^k)\b{A}^{-1} \
& = & (\b{I} - \b{A^{-1}} \b{u} \sum_{k=1}^{\infty}(\b{v^t A^{-1} u})^{k-1}\b{v}) \b{A}^{-1} \
& = & (\b{I} - \b{A^{-1}} \b{u} ~ \frac{\b{v^t}}{1 + \b{v^t A^{-1} u}}) \b{A}^{-1} \
& = & \b{A}^{-1} - \frac{\b{A^{-1} u v^t A^{-1}}}{1 + \b{v^t A^{-1} u}}
\end{eqnarray}
$$

方程式を用いた方法

方程式を用いた方法では以下の $ \b{x} $ 求める過程で目的の式を得ることができる。

$$
\begin{eqnarray}
(\b{A} + \b{uv^t}) \b{x} = \b{b}
\end{eqnarray}
$$

両辺に $ \b{v^t} $ を掛けて $ \alpha = \b{v^t x}$ とおくと

$$
\begin{eqnarray}
\b{Ax} + \b{u}\b{v^tx} & = & \b{b} \
\b{Ax} + \b{u}\alpha & = & \b{b} \
\b{v^t x} & = & \b{v^t} (\b{A^{-1}b} - \b{A^{-1} u \alpha}) \
\alpha & = & \frac{\b{v^t A^{-1} b}}{1+\b{v^t A^{-1} u}}
\end{eqnarray}
$$

これより上記の $ \alpha $ を元の式に代入すると求めたい式が得られる。

$$
\begin{eqnarray}
\b{x} & = & \b{A^{-1}b} - \b{A^{-1} u}\alpha \
\b{x} & = & (\b{A^{-1}} - \frac{\b{A^{-1} u v^t A^{-1}}}{1+\b{v^t A^{-1} u}}) \b{b}
\end{eqnarray}
$$

これより、直接逆行列を掛けた場合と比較して以下が成り立つことが分かる。

$$
(\b{A} + \b{uv^t})^{-1} = \b{A^{-1}} - \frac{\b{A^{-1} u v^t A^{-1}}}{1+\b{v^t A^{-1} u}}
$$

Woodbury の公式

今回は使わないが以下のように Woodbury の公式もノイマン級数を用いた方法で簡単に導出できる。Woodbury の公式は以下の通り。

$$
\begin{eqnarray}
(\b{A} + \b{BCD})^{-1} = \b{A}^{-1} - \b{A^{-1} B} (\b{C^{-1}} + \b{D A^{-1} B})^{-1} \b{D A^{-1}}
\end{eqnarray}
$$

この公式は

• $ \b{A} $ が $ n \times n $
• $ \b{B} $ が $ n \times k $
• $ \b{C} $ が $ k \times k $
• $ \b{D} $ が $ k \times n $

の時に $ k << n $ である場合で $ \b{A}^{-1}, \b{C}^{-1} $ が求められている時、左辺より右辺を計算した方が遥かに計算量が小さくなるため、そのような場面で使われる。

さて、Sherman–Morrison のノイマン級数で導いた時と同様にこの式を証明してみよう。

$$
\begin{eqnarray}
(\b{A} + \b{B C D})^{-1} & = & \{ \b{A} (\b{I} + \b{A^{-1} B C D}) \}^{-1} \
& = & (\b{I} + \b{A^{-1} B C D})^{-1} \b{A^{-1}} \
& = & (\sum_{k=0}^{\infty} (\b{- A^{-1} B C D})^{k}) \b{A^{-1}} \
& = & (\b{I} - \b{A^{-1} B} \{ \sum_{k=1}^{\infty} (\b{- C D A^{-1} B})^{k-1} \} \b{C D} ) ~ \b{A^{-1}} \
& = & \b{A^{-1}} - \b{A^{-1} B}(\b{I} + \b{C D A^{-1} B})^{-1} \b{C D A^{-1}} \
& = & \b{A^{-1}} - \b{A^{-1} B}(\b{C^{-1}} + \b{D A^{-1} B})^{-1} \b{D A^{-1}}
\end{eqnarray}
$$

とこのようにノイマン級数を用いれば比較的簡単に式を導くことが出来る。

Appendix D: ベルヌーイ数の求め方

ディガンマ関数やその微分の漸近展開にはベルヌーイ数が必要である。ベルヌーイ数を求める閉じた形の数式は存在しないので、漸化式を利用して再帰的に値を求める。ベルヌーイ数の母関数の形を思い出すと

$$
\begin{eqnarray}
\frac{x}{e^x-1} & = & \sum_{n=0}^{\infty} \frac{B_n}{n!} x^n \
\end{eqnarray}
$$

より、両辺に $ e^x - 1 $ を掛けることで

$$
\begin{eqnarray}
x & = & \sum_{n=0}^{\infty} \frac{B_n}{n!} (e^x - 1) x^n \
& = & \sum_{n=0}^{\infty} \sum_{k=0}^{\infty} \frac{B_n}{(k+1)! n!} x^{n + k + 1} \
& = & \sum_{l=0}^{\infty} \sum_{n+k=l \substack n, k \ge 0} \frac{B_n}{(l - n + 1)!~n!} x^{l+1} \
& = & \sum_{l=0}^{\infty} (\sum_{n=0}^{l} \frac{1}{(l+1)!} \dbinom{l+1}{n} ~ B_n) x^{l+1} \
\end{eqnarray}
$$

両辺の係数を比較することで、 $ B_0 = 1 $ かつ $ n \ge 2 $ において

$$
\begin{eqnarray}
\sum_{n=0}^{l} \dbinom{l+1}{n} B_n = 0 \
B_l = - \frac{1}{l+1} \sum_{n=0}^{l-1} \dbinom{l+1}{n} B_n
\end{eqnarray}
$$

が成り立つことが分かる。
再帰的に上記の式を計算することで、任意の $ B_n $ の値が分かる。

参考までに ruby で書いたコードと40までのベルヌーイ数の具体的な値を表として残しておく。

def fac(n)
  return 1 if n == 0
  return (1..n).reduce(&:*)
end
def comb(n, k) fac(n) / (fac(n - k) * fac(k)) end
def bernoulli_number(m)
  @memo ||= {}
  return @memo[m] if @memo[m]
  return 1 if m == 0
  return @memo[m] = (0..(m-1)).map {|n|
    comb(m + 1, n) * bernoulli_number(n)
  }.reduce(&:+) * Rational(-1, m + 1)
end
(0..100).each {|x| puts "#{x}: #{bernoulli_number(x)}"}

index	value
0	1
1	-1/2
2	1/6
3	0/1
4	-1/30
5	0/1
6	1/42
7	0/1
8	-1/30
9	0/1
10	5/66
11	0/1
12	-691/2730
13	0/1
14	7/6
15	0/1
16	-3617/510
17	0/1
18	43867/798
19	0/1
20	-174611/330
21	0/1
22	854513/138
23	0/1
24	-236364091/2730
25	0/1
26	8553103/6
27	0/1
28	-23749461029/870
29	0/1
30	8615841276005/14322
31	0/1
32	-7709321041217/510
33	0/1
34	2577687858367/6
35	0/1
36	-26315271553053477373/1919190
37	0/1
38	2929993913841559/6
39	0/1
40	-261082718496449122051/13530