条件付き確率場(CRF)概説 CRF編 - 気まぐれエンジニア日記

前回の記事でCRF以前の多クラス線形識別器まで話をしました. 今回はその続きということで, 本題であるCRFについてみていくことにします.
その前回の記事というのがこれです
kumachan-math.hatenablog.jp
これを前提に話を進めるのでまだの方は先にこっちから読むことをお勧めします.
前置きが長くても仕方ないので早速本題に行きましょう.

1. CRFは何をするか

多クラス線形識別においては, 1つの入力 $x$ に対して, 1つの識別結果 $y$ を与えるのでした. 具体的には, ソフトマックス関数によって定義される"条件付き確率"のようなもの
$\displaystyle p(y|x) := \frac{\exp(\textbf{v}^T\textbf{f}(x, y))}{\displaystyle \sum_{y' \in Y} \exp(\textbf{v}^T\textbf{f}(x, y'))}$
を定義し, 正解ラベル $y^*$ に対するこの値の積による尤度, あるいはその対数を取った対数尤度
$\displaystyle L(Y|X) := \prod_{x \in X} p(y^*|x) = \prod_{x \in X} \frac{\exp(\textbf{v}^T\textbf{f}(x, y^*))}{\displaystyle \sum_{y' \in Y} \exp(\textbf{v}^T\textbf{f}(x, y'))}$
$\displaystyle \log L(Y|X) = \log\left(\prod_{x \in X} p(y^*|x)\right) = \sum_{x \in X} \log p(y^*|x) = \sum_{x \in X} \log\left(\frac{\exp(\textbf{v}^T\textbf{f}(x, y^*))}{\displaystyle \sum_{y' \in Y} \exp(\textbf{v}^T\textbf{f}(x, y'))}\right)$

の最大化問題として定式化されるのでした. けれども, 例えば自然言語や音声データを対象として何らかの識別問題を解きたいと思った時には, 入力 $x$ の識別結果が $x$ の情報のみによって決まると考えるのは安易です. 自然言語であれば単語単体ではなく文脈を考慮したような識別ができた方がいいですし, 音声でも前後の発話からの抑揚の変化だったり速度の変化だったりを考慮できた方が良いですよね.
CRFは, 何らか順序を持ったデータをそれ単体で扱うのではなく系列として扱うということをします. 入力は今まで $x$ だったものが順序を持った複数の入力の系列 $\textbf{x} = x_1x_2, \cdots, x_M$ になるし, 出力も $\textbf{y} = y_1y_2, \cdots, y_M$ になります.
問題は複雑になりましたが, 式の上ではあまり問題はあまり変わっていません. 最大化するべき対数尤度は次のようになります.
$\displaystyle \log L(\textbf{Y}|\textbf{X}) = \log\left(\prod_{\textbf{x} \in \textbf{X}} p(\textbf{y}^*|\textbf{x})\right) = \sum_{\textbf{x} \in \textbf{X}} \log p(\textbf{y}^*|\textbf{x}) = \sum_{\textbf{x} \in \textbf{X}} \log\left(\frac{\exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}^*))}{\displaystyle \sum_{\textbf{y}' \in \textbf{Y}} \exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}'))}\right)$

ここで, 出力に関しては $\textbf{Y} = Y^M$ ですが, 入力は順序関係の束縛があるため真部分集合 $\textbf{X} \subset X^M$ です.
以上です. CRFってつまりはこれだけです. 数式の上では. 大事なことなのでもう一度言います.
数式の上では
これで終わりです. つまり, まだまだ実際に動くシステムには遠いということです. 問題点を挙げましょう.

$\textbf{f}(\textbf{x}, \textbf{y}^*)$ ってさらっと書いてありますが, これはどうやって見つけるのでしょうか. 実験的に頑張るのでもいいかもしれませんが, もう少し理論に裏打ちされたセンスのある見つけ方というものはないのでしょうか
式を複雑にするのは結構ですが, コンピュータに計算させるとなると計算時間の問題は常について回ります. どうにかして高速に計算できないでしょうか.

例えば, 尤度の式の中に出てくる
$\displaystyle Z_\textbf{x} = \sum_{\textbf{y}' \in \textbf{Y}} \exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}'))$
とかは計算が大変そうです. 他にもパラメタ更新を勾配法で行うとき, その更新式の計算は一筋縄では行きそうにありません.

これらの点を順に見ていきましょう.

2. マルコフ性の仮定の導入とCRFの識別

今, 推定したいラベル列を $\textbf{y} = y_1\cdots y_{m-1}y_my_{m+1}\cdots y_M$ とし, これについて各 $y_m$ は $y_{m-1}, x_m$ のみに依存して決まるという仮定をおきましょう. この仮定は問題によって変更をするところですが, ひとまず説明のためにこうしておきます. この仮定の下では次の3種類の特徴を得ることができます.

$y_m$ と $y_{m-1}$ と $x_m$ が全て関係する特徴
$y_m$ と $x_m$ のみが関係する特徴
$y_m$ と $y_{m-1}$ のみが関係する特徴

$y_{m-1}$ の値はすでに決まっていますから, $x_m$ と $y_{m-1}$ のみが関係する特徴は考える必要がありません. これらの特徴をいくつか定義することができれば, それらを各成分にもつベクトルを $\textbf{x}$ のうちの位置 $m$ における局所的な特徴ベクトルと考えることができます. この局所的な特徴ベクトルを $\pi(y_m, y_{m-1}, x_m)$ と書くことにします. $\textbf{x}$ 全体の特徴ベクトルを得るときには, 最も単純にはそれらを足し合わせれば良いですから.
$\displaystyle \textbf{f}(\textbf{x}, \textbf{y}) := \sum_{m = 1}^M \pi(y_m, y_{m-1}, x_m)$
と定義すれば良さそうです.
まずは, 学習のステップは考えないで, 学習が終わった後の識別のステップを考えましょう. $\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y})$ の最大化に注目します. 最大化における変数は他でもなく識別結果のラベル列 $\textbf{y} = y_1\cdots y_{m-1}y_my_{m+1}\cdots y_M$ です.
$\displaystyle \textbf{v}^T \textbf{f}(\textbf{x}, \textbf{y}) = \textbf{v}^T \sum_{m = 1}^M \pi(y_m, y_{m-1}, x_m) = \sum_{m = 1}^M \textbf{v}^T \pi(y_m, y_{m-1}, x_m)$
ですから,
$\displaystyle r_\textbf{x}(m, y_m, y_{m-1}) := \textbf{v}^T \pi(y_m, y_{m-1}, x_m)$
と定義すれば,
$\displaystyle \textbf{v}^T \textbf{f}(\textbf{x}, \textbf{y}) = \sum_{m = 1}^M r_\textbf{x}(m, y_m, y_{m-1})$
とかけます. 改めて問題を明示的に式で書けば,
$\displaystyle \textbf{y}^* = \arg\max_{\textbf{y} \in \textbf{Y}} \sum_{m = 1}^M r_\textbf{x}(m, y_m, y_{m-1})$
です.

これ, どうやってときましょうか. $\textbf{y} \in \textbf{Y} = Y^M$ は組合せ爆発の問題を抱えていて, 時間軸上順番を持った系列を考えており, 各 $y_m$ は $Y$ の元それぞれを選択肢として持っていて, $y_1\cdots y_{m-1}$ でどんなラベル列を選んできたかが $y_m$ の選び方を束縛しないわけです. こうきたら, 典型的な動的計画法...ですよね. DPと分かれば漸化式を書いてみましょう. 時間的な流れ(今まで変数を $m$ と書きました)を一方に, それぞれの選択肢たるラベル $Y$ の元(今まで変数を $y$ と書きました)を他方の次元に持つようなDP行列 $R(m, y)$ を考えれば良さそうです.
$\displaystyle R(1, y) = r_\textbf{x}(1, y, \mathrm{BOS}), ~~~ R(m, y) = \max_{y' \in Y} \{R(m-1, y') +r_\textbf{x}(0, y, y')\}$
ここで, $\mathrm{BOS}$ はBeginning Of Sequenceのことで, 系列の開始を表すダミーのラベルです. あとはこのDP行列を計算して, 最大成分のところからバックトラックすれば, 識別結果のラベル列 $\textbf{y} = y_1\cdots y_M$ が得られます. バックトラック...説明した方が良ければコメントでください.

3. 分配関数の計算

学習のステップにおいて最大化したい尤度の中に
$\displaystyle Z_\textbf{x} = \sum_{\textbf{y}' \in \textbf{Y}} \exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}'))$
という式が出てきました. この関数のことは分配関数と呼んでいます. これを愚直に計算しようと思うと $|Y|^M$ 回もの内積計算が必要になります. つまりラベルの種類の数に関して指数オーダーです. 誰かさんの言葉を借りれば, 恥ずかしいほど遅いを通り越して論外です. いいですか, 論 of the 外です. もう少し効率の良い計算方法はないものだろうか, そう考えるのは自然なことです.
せっかく上で $r_\textbf{x}(m, y_m, y_{m-1})$ を定義したので, それを用いて分配関数の定義を書き直してみましょう.
$\displaystyle Z_\textbf{x} = \sum_{\textbf{y}' \in \textbf{Y}} \exp(\sum_{m = 1}^M r_\textbf{x}(m, y'_m, y'_{m-1})) = \sum_{\textbf{y}' \in \textbf{Y}} \prod_{m = 1}^M \exp(r_\textbf{x}(m, y'_m, y'_{m-1}))$
式の見た目が煩雑になるのが嫌なので,
$\displaystyle \hat{r}_\textbf{x}(m, y'_m, y'_{m-1}) := \exp(r_\textbf{x}(m, y'_m, y'_{m-1}))$
と置きます. この(因数分解と展開の意味で)展開された式を因数分解してみようと思うと
$\displaystyle Z_\textbf{x} = \sum_{\textbf{y'} \in Y^M}^K \prod_{m = 1}^M \hat{r}_\textbf{x}(m, y'_m, y'_{m-1}) = \sum_{y'_1 \in Y} \hat{r}_\textbf{x}(1, y'_1, \mathrm{BOS}) \sum_{y'_2 \in Y} \hat{r}_\textbf{x}(2, y'_2, y'_1) \sum_{y'_3 \in Y} \cdots \sum_{y'_M \in Y} \hat{r}_\textbf{x}(M, y'_M, y'_{M-1})$
となります. 少々煩雑ですが, ただの因数分解なので, ゆっくり式を追えば必ず理解できます(頑張って！). ここまでくれば, この計算を再帰的に行うことができると言うことに気がつくと思います.
$\displaystyle \zeta_\textbf{x}(m, y'_m) = \sum_{y'_1 \in Y} \hat{r}_\textbf{x}(1, y'_1, \mathrm{EOS}) \sum_{y'_2 \in Y} \cdots \sum_{y'_{m-1} \in Y} \hat{r}_\textbf{x}(m-1, y'_{m-1}, y'_{m-2}) \hat{r}_\textbf{x}(m, y'_m, y'_{m-1})$
とおくと, これらは漸化式
$\displaystyle \zeta_\textbf{x}(1, y'_1) = \hat{r}_\textbf{x}(1, y'_1, \mathrm{EOS}), ~~ \zeta_\textbf{x}(m, y'_m) = \sum_{y'_{m-1}\in Y} \zeta_\textbf{x}(m-1, y'_{m-1})\hat{r}_\textbf{x}(m, y'_m, y'_{m-1})$
を満たします. あとはこの漸化式に従って行列を埋めていけば,
$\displaystyle Z_\textbf{x} = \sum_{y'_M \in Y} \zeta_\textbf{x}(M, y'_M)$
として分配関数を求めることができました. 組み合わせ爆発で困ったら, とりあえずDPを検討すればいいみたいですね.

次の章のために, もう一種類の漸化式を考えておきましょう.
$\displaystyle \eta_\textbf{x}(m, y'_{m-1}) :=\sum_{y'_m\in Y} \hat{r}_\textbf{x}(m, y'_m, y'_{m-1})\sum_{y'_{m+1}\in Y} \cdots \sum_{y'_M\in Y} \hat{r}_\textbf{x}(M, y'_M, y'_{M-1})$
とおくと, これらは漸化式
$\displaystyle \eta_\textbf{x}(M, y'_{M-1}) = \sum_{y'_M\in Y} \hat{r}_\textbf{x}(M, y'_M, y'_{M-1}), ~~ \eta_\textbf{x}(m, y'_{m-1}) = \sum_{y'_m\in Y} \hat{r}_\textbf{x}(m, y'_m, y'_{m-1})\eta_\textbf{x}(m+1, y'_m)$
を満たして,
$\displaystyle Z_\textbf{x} = \eta_\textbf{x}(1, \mathrm{EOS})$
となります.

4. CRFのオンライン学習

最尤推定によって学習を進めるにあたって, 謎の変数がまだ1つ残っています. というより, この謎の変数が学習によって調整をするものな訳ですけれどね. そう,今まで定ベクトルとか言ってだましだまし考えていた $\textbf{v}$ です. この $\textbf{v}$ をうまく決めましょう, ただし評価関数は対数尤度
$\displaystyle \log L(\textbf{Y}|\textbf{X}) = \sum_{\textbf{x} \in \textbf{X}}\log p(\textbf{y}^*|\textbf{x}) = \sum_{\textbf{x} \in \textbf{X}} \log\left(\frac{\exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}^*))}{Z_\textbf{x}(\textbf{v})}\right)$
にしましょうというのがCRFの学習そのものです. オンライン学習で勾配法によって $\textbf{v}$ を更新するときの更新式は $\alpha_t$ を学習率として,
$\textbf{v} \leftarrow \textbf{v} - \alpha_t (\nabla l)( \textbf{v})$
とかけます. ここに, $\textbf{x}, \textbf{y}^*$ を新しく追加したデータとして
$\displaystyle l(\textbf{v} | \textbf{x}, \textbf{y}^*) := \log p(\textbf{y}^*|\textbf{x}) = \log\left(\frac{\exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}^*))}{Z_\textbf{x}(\textbf{v})}\right) = \textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}^*) - \log Z_\textbf{x}(\textbf{v})$
です. このとき,
$\displaystyle (\nabla l)( \textbf{v}) = \textbf{f}(\textbf{x}, \textbf{y}^*) - \frac{1}{Z_\textbf{x}(\textbf{v})}(\nabla Z_\textbf{x})( \textbf{v})$
となりますが,
$\displaystyle (\nabla Z_\textbf{x})( \textbf{v}) = \sum_{\textbf{y}' \in \textbf{Y}} \exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}'))\textbf{f}(\textbf{x}, \textbf{y}')$
ですから,
$\displaystyle (\nabla l)( \textbf{v}) = \textbf{f}(\textbf{x}, \textbf{y}^*) - \frac{1}{Z_\textbf{x}(\textbf{v})}\sum_{\textbf{y}' \in \textbf{Y}} \exp(\textbf{v}^T\textbf{f}(\textbf{x}, \textbf{y}'))\textbf{f}(\textbf{x}, \textbf{y}') = \textbf{f}(\textbf{x}, \textbf{y}^*) - \sum_{\textbf{y}' \in \textbf{Y}} p(\textbf{y}'|\textbf{x}) \textbf{f}(\textbf{x}, \textbf{y}')$
と綺麗な形に整理されました.

今から, 第2項を $\zeta_\textbf{x}(m, y'_m)$ と $\eta_\textbf{x}(m, y'_{m-1})$ を使って表現することを考えましょう. もしこれがうまくいけば, 分配関数の計算過程で出てきた副産物を使いまわすことができるので計算時間はグッと短縮されます. 第2項の総和では $Y^M$ の元全てについて考えるのですから, 局所的にある連続した2つのラベルを見たときには, $Y^2$ の全てについて考えていると言えます. 入力として $\textbf{x}$ が与えられたとき位置 $m$ において, 直前のラベルが $y'$ であって位置 $m$ のラベルが $y$ となるソフトマックス関数の意味での確率を $p(y, y'|\textbf{x}, m)$ と書くことにしましょう. すると, 上の第二項は $\pi(y_m, y_{m-1}, x_m)$ を用いて,
$\displaystyle \sum_{\textbf{y}' \in \textbf{Y}} p(\textbf{y}'|\textbf{x}) \textbf{f}(\textbf{x}, \textbf{y}') = \sum_{(y, y') \in Y^2} p(y, y'|\textbf{x}, m) \pi(y, y', x_m)$
とかけます. $p(y, y'|\textbf{x}, m)$ は, 分散関数の値を分母に持って, その計算過程において, (m-1, y')と(m, y)を通る値の総和を分子に持つ比によって書けるので, ある位置 $m$ より先の値の総和を計算している $\eta_\textbf{x}(m+1, y_{m-1})$ とその直前の位置 $m-1$ までの値の総和を計算している $\zeta_\textbf{x}(m-1, y')$ 用いれば,
$\displaystyle p(y, y'|\textbf{x}, m) = \frac{\eta_\textbf{x}(m+1, y)\hat{r}_\textbf{x}(m, y, y')\zeta_\textbf{x}(m-1, y')}{Z_\textbf{x}}$
となります. 晴れて目標達成です.

5. まとめ

CRF概説はここまでとなります. 色々数式をいじりましたが, お気持ちをまとめると次の二点です.

CRFは, 何らか順序を持ったデータをそれ単体で扱うのではなく系列として扱い, 入力系列 $\textbf{x}$ から識別ラベルの系列 $\textbf{y}$ を推定するための手法です.
ソフトマックス関数によって確率のようなものを定義して, それによる尤度の最大化問題として学習ステップを行い, 内積によるスコアの最大化問題として識別ステップをおこないます.

はー長かったですね. お疲れした〜！