$$
\newcommand{\b}[1]{\mathbf{#1}}
\newcommand{\t}[1]{\tilde{#1}}
\newcommand{\max}{\mathop{\rm max}\limits}
\newcommand{\argmax}{\mathop{\rm argmax}\limits}
$$
自然言語処理における条件付き確率場(CRF)の導出
条件付き確率場は単語分割や品詞推定など自然言語処理の様々な場面で利用される手法である。
自分でスクラッチからCRFを実装するために、どのような統計モデルに基いてどのようにして導出するのかをまとめてみた。
問題設定
条件付き確率場は系列 $ x_1x_2...x_N $ が与えられた時にその系列に対するタグや品詞 $ y_1y_2...y_N $ を推測することが出来る。
ここで単語の集合を $ Word $, 品詞の集合を $ WordClass $ とすると、$
x_i \in Word \ \ (\forall i \in {1, 2, ... N}) $ かつ $
y_i \in WordClass \ \ (\forall i \in {1, 2, ... N}) $ である。すなわち 各 $ x_i $ は単語集合の中から選ばれ、 $ y_i $ は品詞集合の中から選ばれる。
また、一つの文における品詞の組合せ方は文長 N の冪乗、 $ |WordClass|^N $ だけ存在する。したがって総当りでこれらの品詞を推定するのは実質的に不可能である。ここで様々な工夫が必要になる。特に動的計画法を上手く活用することがポイントになる。
問題設定(2)（アノテーションの場合）
文章の中から特定の概念を表す文字列（例えば映画の題名や人の人名など）を取得することも同様にCRFを用いると実現できる。この場合、$ x_1x_2...x_N $ は品詞推定の場合と同様に文章である。一方で $ y_1y_2...y_N $ については文章が目的文字列の中にあるか、外にあるかのラベルになる。具体的には、目的文字列の中の文字の場合 "C"、目的文字列の始まりの場合、"S"目的文字列の終了の場合"E"、目的文字列以外の場合"N" のそれぞれの 4 種類の中からラベリングされる。
S で始まり C で続いて E で終了する部分の文字列を抜き出せば、人名や目的の文字列を取り出すことができる。
モデル化
系列 $ \b{x} = {x_1, x_2, ... x_N} $ を観測したという前提のもと品詞列 $ \b{y} = {y_1, y_2, ... y_N} $ を適切にモデル化したい。
そのような品詞列 $\b{y}$ は $ \b{x} $ を与えられた時にエントロピー $ \sum_{\b{x}} P(\b{x} | \b{y}) log P(\b{x} | \b{y}) $ を最大化させることで得られる。このエントロピーを最大化させる確率分布は以下の式で表すことが出来る（この式の導出は今回行わない。あくまで以下の確率分布で表せるという前提で適切な$\b{y}$を求める）
$$
\begin{eqnarray}
P(\b{y} | \b{x}) = \frac{1}{Z_{\b{x}, \b{w}}} exp(\b{w} \cdot \phi(\b{x}, \b{y}))
\end{eqnarray}
$$
$ Z_{\b{x}, \b{w}} $ は正規化のための分配関数で規格化条件を満たすような値に設定する。
$$
\begin{eqnarray}
\sum_{\b{y}} P(\b{y} | \b{x}) = 1 \\
\sum_{\b{y}} \frac{1}{Z_{\b{x}, \b{w}}} exp(\b{w} \cdot \phi(\b{x}, \b{y})) = 1 \\
Z_{\b{x}, \b{w}} = \sum_{\b{y}} exp(\b{w} \cdot \phi(\b{x}, \b{y}))
\end{eqnarray}
$$
また $ \phi(\b{x}, \b{y}) $ は系列 $ \b{x}, \b{y} $ が与えられた時の特徴ベクトルを表す。（$ \phi(\b{x}, \b{y}) $ は小文字で書かれているがベクトルなので注意すること）
さらに仮定として $ x_t, y_t, y_{t-1} (t = 1, 2, ... N) $ は独立ではないのが、その他の変数はそれぞれ独立であると仮定すると
$$
\begin{eqnarray}
\phi(\b{x}, \b{y}) = \sum_{t = 1}^N \phi(x_t, y_t, y_{t-1})
\end{eqnarray}
$$
として表せる。
特徴量
自然言語処理などの品詞推定では $ x_t, y_t, y_{t-1} $ を全て独立とするのではなく $ x_t, y_t $ と $ y_t, y_{t-1} $ をそれぞれ独立であると仮定して特徴量を作る。
$$
\begin{eqnarray}
\eta(x_t, y_t) & = & { \delta(X = x_t)\delta(Y = y_t) \ | \ x_t \in Word, y_t \in WordClass } \\
\varphi(y_t, y_{t-1}) & = & { \delta(Y = y_t)\delta(Y = y_{t-1}) \ | \ y_t, y_{t-1} \in WordClass }
\end{eqnarray}
$$
とすると特徴量ベクトルは次のようになる。
$$
\begin{eqnarray}
\phi(x_t, y_t, y_{t-1}) & = & \eta(x_t, y_t) \cup \varphi(y_t, y_{t-1})
\end{eqnarray}
$$
推論
最尤推定法により、事例集合 $ DataSet = \{ (\b{x}^{(1)}, \b{y}^{(1)}), (\b{x}^{(2)}, \b{y}^{(2)}), ..., (\b{x}^{(D)}, \b{y}^{(D)}) \} $ が与えられた時の確率 $ \prod_{i=1}^{D} P(\b{x^{(i)}} | \b{y^{(i)}}) $ が最大になるようにパラメーター $ \b{w} $ を定める。
すなわち求めべきパラメーター $ \b{w^{(op)}} $ は
$$
\begin{eqnarray}
\b{w^{(op)}} &=& \argmax_{\b{w}} \prod_{i=1}^{D} P(\b{x^{(i)}} | \b{y^{(i)}}) \\
&=& \argmax_{\b{w}} \ log (\prod_{i=1}^{D} P(\b{x^{(i)}} | \b{y^{(i)}})) \\
&=& \argmax_{\b{w}}\ L(\b{w}) \\
\end{eqnarray}
$$
である。このようにしてパラメーター $ \b{w^{(op)}} $ 求めることが出来る。
ここで $ L(\b{w}) = log(\prod_{i=1}^{D} P(\b{x^{(i)}} | \b{y^{(i)}})) $ は目的関数として新しく定義した。
さらにこの式は正則化項を含めて次のように具体的に表せる。
$$
\begin{eqnarray}
L(\b{w}) &=& log \prod_{i=1}^D P(\b{y^{(i)}} | \b{x^{(i)}}) - \frac{C}{2} ||\b{w}||^2 \\
&=& \sum_{i=1}^D ( log { \frac{1}{Z_{\b{x}^{(i)}, \b{w}}} exp(\b{w} \cdot \sum_{t=1}^N \phi(x_t^{(i)}, y_{t}^{(i)}, y_{t-1}^{(i)})) } ) - \frac{C}{2} ||\b{w}||^2 \\
&=& (\sum_{i=1}^D (\sum_{t=1}^N \b{w} \cdot \phi(x_t^{(i)}, y_{t}^{(i)}, y_{t-1}^{(i)})) - log Z_{\b{x}^{(i)}, \b{w}}) - \frac{C}{2} ||\b{w}||^2
\end{eqnarray}
$$
$ \b{w^{(op)}} $ をどのように求めるかについては後ほど議論することにして今回は $ \b{w^{(op)}} $ を仮に推定できた場合どのようにして $ \b{y} $ を求めるかについて説明する。
新しく文章 $ \b{x} \ (x_1x_2...x_N) $ が与えられた際の品詞を特定したい。
これには条件付き確率を最大にするような品詞系列 $ \b{y^{(op)}} $ を求めればよい。数式で表すと、$ \b{y^{(op)}} = \argmax_{\b{y}} P(\b{y} | \b{x}) $ である。しかしながら最初に述べたようにこれを直接求めるのは $ |WordCount|^N $ の場合だけ組合せが必要なので直接求めることは難しい。そこで動的計画法の一つである viterbi アルゴリズムを用いて計算する必要がある。
$$
\begin{eqnarray}
log({P(\b{y} | \b{x})}) & = \sum_{t=1}^N \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) + Const
\end{eqnarray}
$$
であるため、動的計画法を用いて次のように効率的に計算できる。
$$
\begin{eqnarray}
f(t, y_t) & = & \max_{y_1, y_2, ... y_{t-1}} \ log({P(\b{y} | \b{x})}) \\
& = & \max_{y_1, y_2, ... y_{t-1}} \sum_{t'=1}^t \b{w^{(op)}} \cdot \phi(x_{t'}, y_{t'}, y_{t'-1}) \\
& = & \max_{y_{t-1}} \ \max_{y_1, y_2, ... y_{t-2}} \sum_{t'=1}^{t-1} \b{w^{(op)}} \cdot \phi(x_{t'}, y_{t'}, y_{t'-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) \\
& = & \max_{y_{t-1}} \ f(t - 1, y_{t-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) \\
\end{eqnarray}
$$
これらをまとめると次のようになる。
$$
\begin{eqnarray}
f(t, y_t) & = & \max_{y_{t-1}} \ f(t - 1, y_{t-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) \\
g(t, y_t) & = & \argmax_{y_{t-1}} \ g(t - 1, y_{t-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1})
\end{eqnarray}
$$
$$ t = 1, 2, ... N $$
$ y_1^{(op)}, y_2^{(op)}, ... y_N^{(op)} $ の値については $ g(t, y_t) \ { t = N-1, ... 2, 1} $ をそれぞれあとから辿ればよい。
$$
\begin{eqnarray}
y_N^{(op)} &=& \argmax_{y_N} g(N, y_N) \\
y_t^{(op)} &=& g(t + 1, y_{t + 1}^{(op)}) \ \ { t = N-1, ... 2, 1} \\
\end{eqnarray}
$$
パラメーターの推定と勾配の計算
目的関数 $ L(\b{w}) $ の最大値を計算するために最急勾配法を用いて更新を行う。
$$
\begin{eqnarray}
\b{w^{(t)}} = \b{w^{(t-1)}} + \epsilon \nabla_{\b{w}} L(\b{w^{(t-1)}})
\end{eqnarray}
$$
具体的に $ L(\b{w}) $ の微分を計算すると
$$
\begin{eqnarray}
\nabla_{\b{w}} L(\b{w}) = \sum_{i=1}^D (\sum_{t=1}^N \phi(x_t^{(i)}, y_{t}^{(i)}, y_{t-1}^{(i)}) - \frac{1}{Z_{\b{w}, \b{x^{(i)}}}} \nabla_{\b{w}} Z_{\b{w}, \b{x^{(i)}}}) - C \b{w}
\end{eqnarray}
$$
分配関数の微分は $ \nabla_{\b{w}} Z_{\b{w}, \b{x^{(i)}}} $ は規格化条件
$$
\begin{eqnarray}
Z_{\b{w}, \b{x}} = \sum_{\b{y}} P(\b{y} | \b{x}) = \sum_{y_1, y_2, ... y_N} exp(\sum_{t=1}^N \b{w} \cdot \phi(x_t, y_t, y_{t-1}))
\end{eqnarray}
$$
より
$$
\begin{eqnarray}
\frac{\nabla_{\b{w}} Z_{\b{w}, \b{x^{(i)}}}}{Z_{\b{w}, \b{x^{(i)}}}} & = & \frac{1}{Z_{\b{w}, \b{x^{(i)}}}} \sum_{y_1, y_2, ... y_N} \nabla_{\b{w}} exp(\sum_{t=1}^N \b{w} \cdot \phi(x_t^{(i)}, y_t, y_{t-1})) \\
& = & \sum_{y_1, y_2, ... y_N} \phi(x_t^{(i)}, y_t, y_{t-1}) \frac{exp(\sum_{t=1}^N \b{w} \cdot \phi(x_t^{(i)}, y_t, y_{t-1}))}{Z_{\b{w}, \b{x^{(i)}}}} \\
& = & \sum_{y_t, y_{t-1}} \phi(x_t^{(i)}, y_t, y_{t-1}) \sum_{y_1, y_2, ... y_{t-2}, y_{t+1}, ... y_{N}} \frac{exp(\sum_{t=1}^N \b{w} \cdot \phi(x_t^{(i)}, y_t, y_{t-1}))}{Z_{\b{w}, \b{x^{(i)}}}} \\
& = & \sum_{y_t, y_{t-1}} \phi(x_t^{(i)}, y_t, y_{t-1}) \ P(y_t, y_{t-1} | x_t^{(i)})
\end{eqnarray}
$$
これより目的関数の勾配は
$$
\begin{eqnarray}
\nabla_{\b{w}} L(\b{w}) = \sum_{i=1}^D \sum_{t=1}^N (\phi(x_t^{(i)}, y_{t}^{(i)}, y_{t-1}^{(i)}) - \sum_{y_t, y_{t-1}} \phi(x_t^{(i)}, y_t, y_{t-1}) \ P(y_t, y_{t-1} | x_t^{(i)})) - C \b{w}
\end{eqnarray}
$$
で表せることが分かる。ここで問題は $ P(y_t, y_{t-1} | x_t^{(i)}) $ をどのようにして計算するかだ。直接全て足しあわせようとすると $ |WordClass|^{(N-2)} $ の計算が必要になる。そこで、ここでも動的計画法を上手く利用して計算回数を減らす。
$ \psi(y_t, y_{t-1}) = exp(\b{w} \cdot \phi(x_t, y_t, y_{t-1})) $$ とおくと $$ P(\b{y} | \b{x}) = \prod_{t=1}^N \frac{1}{Z_{\b{x, w}}} \psi(y_t, y_{t-1}) $ と表せる。しがたって
$$
\begin{eqnarray}
P(y_t, y_{t-1} | x_t) & = & \sum_{y_1, y_2, ... y_{t-2}} \sum_{y_{t+1}, ... y_{N}} \prod_{t'=1}^N \frac{1}{Z_{\b{x, w}}} \psi(y_{t'}, y_{t'-1}) \\
& = & \frac{\psi(y_t, y_{t-1})}{Z_{\b{x, w}}} \sum_{y_1, y_2, ... y_{t-2}} \sum_{y_{t+1}, ... y_{N}} \prod_{\substack{ t=1 \\ t \neq t'}}^N \psi(y_{t'}, y_{t'-1}) \\
& = & \frac{\psi(y_t, y_{t-1})}{Z_{\b{x, w}}} (\sum_{y_1, y_2, ... y_{t-2}} \prod_{t'=1}^{t-1} \psi(y_{t'}, y_{t'-1})) (\sum_{y_{t+1}, ... y_{N}} \prod_{t'=t+1}^N \psi(y_{t'}, y_{t'-1}))\\
& = & \frac{1}{Z_{\b{x, w}}} \psi(y_t, y_{t-1}) \alpha(t - 1, y_{t-1}) \beta(t, y_t)
\end{eqnarray}
$$
ここで
$$
\begin{eqnarray}
\alpha(t - 1, y_{t-1}) & = & \sum_{y_1, y_2, ... y_{t-2}} \prod_{t'=1}^{t-1} \psi(y_{t'}, y_{t'-1}) \\
\beta(t, y_t) & = & \sum_{y_{t+1}, ... y_{N}} \prod_{t'=t+1}^N \psi(y_{t'}, y_{t'-1})
\end{eqnarray}
$$
とそれぞれおいた。 $$ \alpha(t - 1, y_{t-1}), \beta(t+1, y_{t+1}) $$ もそれぞれ再帰的に計算できて
$$
\begin{eqnarray}
\alpha(t - 1, y_{t-1}) & = & \sum_{y_1, y_2, ... y_{t-2}} \prod_{t'=1}^{t-1} \psi(y_{t'}, y_{t'-1}) \\
& = & \sum_{y_{t-2}} \{ \sum_{y_1, y_2, ... y_{t-3}} (\prod_{t'=1}^{t-2} \psi(y_{t'}, y_{t'-1})) \} \ \psi(y_{t-1}, y_{t-2}) \\
& = & \sum_{y_{t-2}} \alpha(t-2, y_{t-2}) \ \psi(y_{t-1}, y_{t-2}) \\
\end{eqnarray}
$$
同様にして $ \beta(t, y_{t}) $ についても再帰的に計算すると
$$
\begin{eqnarray}
\beta(t, y_t) & = & \sum_{y_{t+1}, y_{t+2}, ... y_N} \prod_{t'=t+1}^{N} \psi(y_{t'}, y_{t'-1}) \\
& = & \sum_{y_{t+1}} \{ \sum_{y_{t+2}, y_{t+3}, ... y_N} \prod_{t'=t+2}^{N} \psi(y_{t'}, y_{t'-1}) \} \psi(y_{t+1}, y_t) \\
& = & \sum_{y_{t+1}} \beta(t+1, y_{t+1}) \psi(y_{t+1}, y_t) \\
\end{eqnarray}
$$
となる。
さらに分配関数も次のように $ \alpha(N, y_N) $ や $ \beta(1, y_1) $ から計算できる。
$$
\begin{eqnarray}
Z_{\b{x}, \b{w}} & = & \sum_{\b{y}} P(\b{y} | \b{x}) \\
& = & \sum_{y_1, y_2, ... y_N} \prod_{t=1}^N exp(\b{w} \cdot \psi(y_t, y_{t-1})) \\
& = & \sum_{y_N} \alpha(N, y_N) \\
& = & \sum_{y_1} \beta(1, y_1) \\
\end{eqnarray}
$$
実装用のまとめ
パラメーターの学習方法はまとめると以下のようになる。
$$
\begin{eqnarray}
\b{w^{(t)}} & = & \b{w^{(t-1)}} + \epsilon \nabla_{\b{w}} L(\b{w^{(t-1)}}) \\
\nabla_{\b{w}} L(\b{w}) & = & \sum_{i=1}^D \sum_{t=1}^N (\phi(x_t^{(i)}, y_{t}^{(i)}, y_{t-1}^{(i)}) - \sum_{y_t, y_{t-1}} \phi(x_t^{(i)}, y_t, y_{t-1}) \ P(y_t, y_{t-1} | x_t^{(i)})) - C \b{w} \\
P(y_t, y_{t-1} | x_t) & = & \frac{1}{Z_{\b{x, w}}} \psi(y_t, y_{t-1}) \alpha(t-1, y_{t-1}) \beta(t, y_t) \\
\psi(y_t, y_{t-1}) & = & exp(\b{w} \cdot \phi(x_t, y_t, y_{t-1})) \\
\alpha(t, y_t) & = & \sum_{y_{t-1}} \alpha(t-1, y_{t-1}) \ \psi(y_t, y_{t-1}) \\
\beta(t, y_t) & = & \sum_{y_{t+1}} \beta(t+1, y_{t+1}) \psi(y_{t+1}, y_t) \\
Z_{\b{x}, \b{w}} & = & \sum_{y_N} \alpha(N, y_N)
\end{eqnarray}
$$
予測は以下のようになる。
$$
\begin{eqnarray}
f(t, y_t) & = & \max_{y_{t-1}} \ f(t - 1, y_{t-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) \ \ { t = 1, 2, ... N } \\
g(t, y_t) & = & \argmax_{y_{t-1}} \ g(t - 1, y_{t-1}) + \b{w^{(op)}} \cdot \phi(x_t, y_t, y_{t-1}) { t = 1, 2, ... N } \\
y_N &=& \argmax_{y_N} g(N, y_N) \\
y_t &=& g(t + 1, y_{t + 1}) \ \ { t = N-1, ... 2, 1} \\
\end{eqnarray}
$$
デバッグ方法
数値解析用のプログラムをデバッグする時は基本的な公式や式の性質が成り立つかを確認してデバッグを行う。
正しくプログラムが動作していれば以下の等式がそれぞれ成り立つはずである。
$$
\begin{eqnarray}
\sum_{y_{t-1}, y_t} P(y_t, y_{t-1} | x_t) & = & 1 \\
L(\b{w + \delta w}) - L(\b{w}) & \simeq & \delta \b{w} \cdot \nabla_{\b{w}} L(\b{w})
\end{eqnarray}
$$
overflow を避けるための計算手法
単純に (48) - (54) を実装すると数百の系列を学習させるだけで計算機上で扱える桁が溢れてしまう。なぜならば(52), (53)の$ \alpha, \beta $がおおよそ系列のべき乗で増していくため、一般的に実数を扱う場合の倍精度（double型）の上限である $ 2^{2^{10}} \simeq 10^{308.3 }\ $ を簡単に超えてしまうからである。しがって、何らかの方法でこれを避けなければならない。
log スケールで計算する方法とスケール因子を利用する方法と 2 通りの方法を紹介する
log スケールで計算する方法
桁が溢れる部分は $ \alpha, \beta$ の計算なのでこの部分だけ log の世界で扱ってあげればよい。
$$
\begin{eqnarray}
\t{\alpha}(t, y_t) & = & log (\alpha(t, y_t)) \\
\t{\beta}(t, y_t) & = & log (\beta(t, y_t))
\end{eqnarray}
$$
このように改めて $ \t{\alpha}, \t{\beta} $ を定義すると、(50)-(54) の漸化式はそれぞれ以下のように表せる。
$$
\begin{eqnarray}
P(y_t, y_{t-1} | x_t) & = & \frac{1}{\sum_{yn}\alpha(N, y_N)} \psi(y_t, y_{t-1}) \alpha(t-1, y_{t-1}) \beta(t, y_t) \\
& = & \psi(y_t, y_{t-1}) \ exp \{ \t{\beta}(t, y_t) + \t{\alpha}(t-1, y_{t-1}) - \t{\alpha}(N, y_N) \} \\
\t{\alpha}(t, y_t) & = & log \{ \sum_{y_{t-1}} \alpha(t - 1, y_{t-1}) \psi(y_t, y_{t-1}) \} \\
& = & log { \sum_{y_{t-1}} exp(\t{\alpha}(t-1, y_{t-1})) \psi(y_t, y_{t-1}) } \\
\t{\beta}(t, y_t) & = & log \{ \sum_{y_{t+1}} \beta(t+1, y_{t+1}) \psi(y_{t+1}, y_t) \} \\
& = & log \{ \sum_{y_{t+1}} exp(\t{\beta}(t+1, y_{t+1})) \psi(y_{t+1}, y_t) \} \\
\end{eqnarray}
$$
ここで漸化式にあらわれている $ log \sum exp $ (一般的に logsumexp と呼ばれている)でも log の中の exp の計算において overflow や underflow が発生する可能性があるが、以下のように計算するとそれを上手く防ぐことができる。
logsumxexp の計算方法
1つ目は最大値で引く方法
$$
\begin{eqnarray}
log \sum_{i=1}^N e^{x_i} & = & log \{ \sum_{i=1, i \ne i'} e^{x_i} + e^{x_{i'}^{(max)}} \} \\
& = & log \{ ~ exp(x_{i'}^{(max)}) (\sum_{i=1} exp(x_i-x_{i'}^{(max)})) \} \\
& = & x_{i'}^{(max)} + log \{ ~ \sum_{i=1} exp(x_i-x_{i'}^{(max)}) ~ \} \\
\end{eqnarray}
$$
$ x_{i'}^{(max)} $ は $ x_1, x_2 ... x_N $ の中で最大値を取る値である。このようにすることで $ x_i - x_{i'}^{(max)} \leq 1.0 $ から指数部が $ e $ を超えることはなくオーバーフローを避けることができる。
2つ目は再帰的に最小値を求める方法
$$
\begin{eqnarray}
g(n) & = & log \{ \sum_{i=1}^n e^{x_i} \} \\
& = & log \{ \sum_{i=1}^{n-1} e^{x_i} + e^{x_n} \} \\
& = & log \{ exp(log(\sum_{i=1}^{n-1} e^{x_i})) + e^{x_n} \} \\
& = & f(g(n-1), x_n) \\
\end{eqnarray}
$$
ただしここで
$$
\begin{eqnarray}
f(x, y) = log \{ e^x + e^y \}
\end{eqnarray}
$$
と定義した。
$ f(x, y) $ は次のように計算すると桁あふれが起こらない。
$$
\begin{eqnarray}
f(x, y) = \begin{cases}
y + log(1 + e^{x-y}) & (x \leq y) \\
x + log(1 + e^{y-x}) & (x > y)
\end{cases}
\end{eqnarray}
$$
このようにすることで exponential の桁の部分が常に 1 以下になるので桁溢れを防ぐことができる。
具体的にこの方法をbackword-forwardアルゴリズムに当てはめると以下の漸化式を得る
$$
\begin{eqnarray}
P(y_t, y_{t-1} | x_t) & = & \psi(y_t, y_{t-1}) \ exp \{ \t{\beta}(t, y_t) + \t{\alpha}(t-1, y_{t-1}) - Z_{N} \} \\
Z_{N} & = & \t{\alpha}(N, y_N^{(max)}) + log \sum_{y_{N}} exp(\t{\alpha}(N, y_N) - \t{\alpha}(N, y_N^{(max)})) \\
\t{\alpha}(t, y_t) & = & \t{\alpha}(t-1, y_{t-1}^{(max)}) + log \{ \sum_{y_{t-1}} exp(\t{\alpha}(t-1, y_{t-1}) - \t{\alpha}(t-1, y_{t-1}^{(max)})) \psi(y_t, y_{t-1}) \} \\
\t{\beta}(t, y_t) & = & \t{\beta}(t+1, y_{t+1}^{(max)}) + log \{ \sum_{y_{t+1}} exp(\t{\beta}(t+1, y_{t+1})-\t{\beta}(t+1, y_{t+1}^{(max)})) \psi(y_{t+1}, y_t) \} \\
\end{eqnarray}
$$
ここで以下の不等式を満たすように $ t' $ を設定する
$$
\begin{eqnarray}
\t{\alpha}(t-1, y_{t-1}^{(max)}) \geq \t{\alpha}(t-1, y_{t-1}), ~~ \forall y_{t-1} \in WordClass \\
\t{\beta}(t+1, y_{t+1}^{(max)}) \geq \t{\beta}(t+1, y_{t+1}), ~~ \forall y_{t+1} \in WordClass
\end{eqnarray}
$$
スケール因子法
logsumexp を用いる方法では本来必要ではない $ exp $ と $ log $ という比較的重い演算が途中で入るためナイーブな実装と比較して非常に遅くなる。したがってできるだけこのような重い演算が入らないようにしたい。
そのための方法の一つがスケール因子法である。
中心的なアイデアは各時刻tにおいて $ \alpha(t, y_t) $ 適切にスケーリングして$ y_t $ に対して規格化することである。すなわち適当にスケールした変数をそれぞれ $ \t{\alpha}(t, y_t), \t{\beta}(t, y_t) $ として
$$
\begin{eqnarray}
\sum_{y_t} \t{\alpha}(t, y_t) = 1 \\
\sum_{y_t} \t{\beta}(t, y_t) = 1
\end{eqnarray}
$$
を満たすように $ \t{\alpha}, \t{\beta} $ をそれぞれ定めることで $ \alpha, \beta $ の桁を必要以上に大きく（または小さく）することを防ぐ。
比例定数を $ c_t^{(\alpha)}, c_t^{(\beta)} $ とそれぞれする。$ \t{\alpha}, \t{\beta} $ どちらもアルゴリズムはほとんど同じなので以下では主に $ \t{\alpha} $ の規格化定数の求め方と式変形について説明する。 $ \alpha $$ は $$ \t{\alpha} $ に比例するので
$$
\begin{eqnarray}
\t{\alpha}(t, y_t) = c_t^{(\alpha)} \alpha(t, y_t)
\end{eqnarray}
$$
また、比例定数は規格化条件より次の性質を満たさなければならない。
$$
\begin{eqnarray}
\sum_{y_t} c_t^{(\alpha)} \alpha(t, y_t) & = & 1 \\
c_t^{(\alpha)} & = & (\sum_{y_t} \alpha(t, y_t))^{-1}
\end{eqnarray}
$$
漸化式は以下のように表せる。
$$
\begin{eqnarray}
\alpha(t, y_t) = \frac{1}{c_{t-1}^{(\alpha)}} \sum_{y_{t-1}} \t{\alpha}(t-1, y_{t-1}) \psi(y_t, y_{t-1}) \\
\end{eqnarray}
$$
これを (80), (82) にそれぞれ代入すると目的の $ \t{\alpha} $ に対する更新式を得る。
$$
\begin{eqnarray}
c_{t}^{(\alpha)} & = & \frac{c_{t-1}^{(\alpha)}}{ \sum_{y_t} \{ \sum_{y_{t-1}} \t{\alpha}(t-1, y_{t-1}) \psi(y_t, y_{t-1}) \} } \\
\t{\alpha}(t, y_t) & = & \frac{c_t^{(\alpha)}}{c_{t-1}^{(\alpha)}} \sum_{y_{t-1}} \t{\alpha}(t-1, y_{t-1}) \psi(y_t, y_{t-1}) \\
\end{eqnarray}
$$
ただし実はこれでも $ c_t^{(\alpha)} $ の式においてアンダーフローが発生する。なぜならば $ \alpha(t, y_t) $ は $ k^t $ というオーダーで大きくなっていくので比例定数もそれに比例してそれを規格化させるために小さくならなければならないからである。したがって改めて次のような変数を定義する。
$$
\begin{eqnarray}
\t{c_t}^{(\alpha)} & = & \frac{c_t^{(\alpha)}}{c_{t-1}^{(\alpha)}} \\
c_0^{(\alpha)} & = & 1 \\
\alpha(t, y_t) & = & (\prod_{t'=1}^t \t{c_{t'}}^{(\alpha)}) \t{\alpha}(t, y_t) \\
\end{eqnarray}
$$
とおくと(84)-(85)の漸化式はそれぞれ次のように表せる。
$$
\begin{eqnarray}
\t{c_t}^{(\alpha)} & = & \frac{c_t^{(\alpha)}}{c_{t-1}^{(\alpha)}} \\
c_0^{(\alpha)} & = & 1 \\
\alpha(t, y_t) & = & (\prod_{t'=1}^t \t{c_{t'}}^{(\alpha)}) \t{\alpha}(t, y_t) \\
\end{eqnarray}
$$
初期値は以下のようにして与える
$$
\begin{eqnarray}
\t{c_{t}}^{(\alpha)} & = & \frac{1}{ \sum_{y_t} { \sum_{y_{t-1}} \t{\alpha}(t-1, y_{t-1}) \psi(y_t, y_{t-1}) } } \\
\t{\alpha}(t, y_t) & = & \t{c_{t}}^{(\alpha)} \sum_{y_{t-1}} \t{\alpha}(t-1, y_{t-1}) \psi(y_t, y_{t-1}) \\
\end{eqnarray}
$$
同様にして $ \t{\beta} $ の更新式も次のようになる
$$
\begin{eqnarray}
\t{c_t}^{(\beta)} & = & \frac{c_{t+1}^{(\beta)}}{c_t^{(\beta)}} \\
\t{c_N}^{(\beta)} & = & c_N^{(\beta)} = \frac{1}{\sum_{y_{N}, y_{N1+1}} \psi(y_{N+1}, y_N)} \\
\t{\beta}(N, y_N) & = & \t{c_N}^{(\beta)} \sum_{y_{N+1}} \psi(y_{N+1}, y_N) \\
\t{c_{t}}^{(\beta)} & = & \frac{1}{ \sum_{y_t} \{ \sum_{y_{t+1}} \t{\beta}(t+1, y_{t+1}) \psi(y_{t+1}, y_t) \} } \\
\t{\beta}(t, y_t) & = & \t{c_t}^{(\beta)} \sum_{y_{t+1}} \t{\beta}(t+1, y_{t+1}) \psi(y_{t+1}, y_t) \\
\end{eqnarray}
$$
周辺分布はこれらの結果を用いて以下のように計算する。
$$
\begin{eqnarray}
P(y_t, y_{t-1} | \b{x}) & = & \frac{1}{\prod_{t'=1}^N \t{c_{t'}}^{(\alpha)}} \phi(y_t, y_{t-1}) \prod_{t'=1}^{t-1} \t{c_{t'}}^{(\alpha)} \t{\alpha}(y_t, y_{t-1})\prod_{t'=t+1}^N \t{c_{t'}^{(\beta)}} \beta(t, y_t) \\
& = & \frac{1}{\t{c_t^{(\alpha)}}} (\prod_{t'=t+1}^N \frac{\t{c_{t'}^{(\beta)}}}{\t{c_{t'}^{(\alpha)}}}) ~ \phi(y_t, y_{t-1}) \t{\alpha}(y_t, y_{t-1}) \t{\beta}(t, y_t)
\end{eqnarray}
$$