続・わかりやすいパターン認識ひとり読書会 2 - 渡邊哲司（Satoshi Watanabe）の日記

前回に引き続きディリクレ過程混合モデルの考察です。今回は具体的にアルゴリズムをどのように組めばよいのかを考えます。テキストの続パタにも十分詳細にアルゴリズムが載っていますがやはり細かい部分では行間が空いていますので、その部分を埋めることが一つの目的です。

続パタ読書会 1:ディリクレ過程混合モデル
続パタ読書会 2:ディリクレ過程混合モデルのアルゴリズム
 続パタ読書会 3:実験
 続パタ読書会 4:実装
 続パタ読書会 5:事前確率P(s)の妥当性 1
続パタ読書会 6:事前確率P(s)の妥当性 2
続パタ読書会 7:ベル数 1
続パタ読書会 8:ベル数 2
続パタ読書会 9:P(s_k|x_k, s_-k, θ)の導出
 続パタ読書会 10:演習問題12.5の計算
 続パタ読書会 11:演習問題12.6の計算 1
続パタ読書会 12:演習問題12.6の計算 2

ディリクレ過程混合モデルのアルゴリズム

概要

ディリクレ過程混合モデルの実装として必要な処理は大まかに

1. 初期化
　(1.) ハイパーパラメータの設定
　(2.) パターンの入力
　(3.) その他のパラメータの設定
　(4.) クラスタ構造の変更に対して不変である式の計算
　(5.) 対数事後確率の初期値の計算
2. 学習
　(1.) $\boldsymbol{s}$ の更新
　(2.) $\boldsymbol{\theta}$ の更新
　(3.) 対数事後確率の計算
3. 結果取得

という感じです。では各々の処理を具体的にどのようにするかを考えていきます。

初期化

1. ハイパーパラメータの設定

学習を開始するために必要なパラメータを事前に設定します。具体的に事前に設定が必要なパラメータは

$\alpha$ ：集中度
$\beta$ ：正規ウィシャート分布のパラメータ
$\nu$ ：自由度
$\mathbf{S}$ ： $d\times d$ 正定値対象行列

の四つです。

2. パターンの入力

続いてパターン集合 $\boldsymbol{x}=\{x_1, \cdots, x_n\}$ を入力します。パターンが入力されることによって同時にパターン総数と特徴空間の次元数も判明しますのでそれらも設定します。つまり

$\boldsymbol{x}$ ：パターン集合
$n$ ：パターン総数
$d$ ：特徴空間の次元数

の3つを設定します。

3. その他のパラメータの設定

パターン集合 $\boldsymbol{x}$ を得たことで設定が可能となる残ったパラメータをすべて初期化します。1.で設定したパラメータはパターンの入力なしに設定が可能なもの、2.で設定したパラメータはパターンが入力されれば同時に確定するパラメータですが、ここで設定するパラメータは事前にどのように設定するかを決めておく必要があるもののうちパターン集合が入力されなければ設定できないパラメータです。ここでもテキストに倣い初期状態では全パターンが同一のクラスタに所属する、つまりクラスタ数は1から始めるものとし、 $\mu_0$ は全パターンの標本平均 $\mu_{all}$ とします。具体的には下記のように設定します。

$\boldsymbol{s}=\{s_1, \ldots, s_n\}$ はすべて $1$

全パターンが同一クラスタのため所属クラスタはすべて $1$ とします。実装では言語(python)上のインデックスが $0$ から始まるためにすべて $0$ としていますが、意味的には全パターンがクラスタ $\omega_1$ に所属することを表していることに変わりはありません。

$c=1$

やはり全パターンが同一クラスタのためクラスタ数 $c$ は $1$ とします。

$n_1=n$

$i$ 番目のクラスタ（つまり $\omega_i$ ）に所属しているパターン数を $n_i$ で表します。ここではクラスタ数が一つしかなく、全パターンがその一つのクラスタに所属しているため、 $n_1$ のみが存在し、所属パターン数はパターン総数である $n$ となります。

$\mu_0=\mu_{all}$

テキストに倣いハイパーパラメータ $\mu_0$ は全パターンの標本平均 $\mu_{all}$ で初期化します。

$\mu_1=\mu_{all}$
$\boldsymbol{\Lambda}_1=\boldsymbol{\Lambda}_{all}$

全パターン同一クラスタのため $\mu_1$ は $\mu_0$ と同じく全パターンの標本平均 $\mu_{all}$ で初期化します。また、精度行列 $\boldsymbol{\Lambda}_1$ は全パターンより算出した分散共分散行列 $\boldsymbol{\Sigma}_{all}$ の逆行列で初期化します。つまり $\boldsymbol{\Lambda}_1=\boldsymbol{\Lambda}_{all}=\boldsymbol{\Sigma}_{all}^{-1}$ です。

4. クラスタ構造の変更に対して不変である式の計算

学習を進めていくとクラスタ構造は変更されますが、クラスタ構造の変更に対して不変なパラメータのみを持つ式の計算内容は変わりません。ですのでここではそれらを事前に計算しておきます。具体的には対数事後確率 $\log{v}$ の項の一つであるイーウェンスの抽出公式

$\begin{eqnarray} \log{P_E(n_1, \ldots , n_c)} &=& \log{\displaystyle\frac{\alpha^c\displaystyle\prod_{i=1}^{c}{(n_i-1)!}}{\text{AF}(\alpha, n)}}\\ &=& c\log{\alpha}+\displaystyle\sum_{i=1}^{c}{\sum_{l=1}^{n_i-1}{\log{l}}}-\log{\text{AF}(\alpha, n)} \end{eqnarray}$

における最初の項と最後の項を計算しておきます。最初の項の $\log{\alpha}$ の部分は集中度 $\alpha$ のみに依存しておりクラスタ構造に対し不変です。また最後の項の $\log{\text{AF}(\alpha, n)}$ は集中度 $\alpha$ とパターン総数 $n$ のみに依存するためやはりクラスタ構造に対し不変です。後者の $\log{\text{AF}(\alpha, n)}$ は実際の計算を書き下ろせば

$\begin{eqnarray} \log{\text{AF}(\alpha, n)} &=& \log{\alpha(\alpha+1)\cdots(\alpha+n-1)}\\ &=& \displaystyle\sum_{k=0}^{n-1}{\log{(\alpha+k)}}\\ \end{eqnarray}$

となりますのでこれを計算しておきます。

5. 対数事後確率の初期値の計算

最後に対数事後確率の初期値を計算します。具体的に書き下ろせば

$\begin{eqnarray} \log{v} &=& \log{\left(P(s)\displaystyle\prod_{i=1}^{c}{\left(p(\theta_i)\displaystyle\prod_{x_k\in\omega_i}{p(x_k\vert\theta_i)}\right)}\right)}\\ &=& \log{P(s)}+\displaystyle\sum_{i=1}^{c}{\log{p(\theta_i)}}+\displaystyle\sum_{i=1}^{c}\sum_{x_k\in\omega_i}\log{p(x_k\vert\theta_i)}\\ &=& \log{P_E(n_1, \ldots , n_c)}\\ && +\displaystyle\sum_{i=1}^{c}{\log{\mathcal{N}(\mu_i;\mu_0, (\beta\mathbf{\Lambda}_i)^{-1})\mathcal{W}(\mathbf{\Lambda}_i;\nu, \mathbf{S})}}\\ && +\displaystyle\sum_{i=1}^{c}\sum_{x_k\in\omega_i}{\log{\mathcal{N}(x_k;\mu_i, \mathbf{\Lambda}_i^{-1})}}\\ &=& c\log{\alpha}+\displaystyle\sum_{i=1}^{c}{\sum_{l=1}^{n_i-1}{\log{l}}}-\log{\text{AF}(\alpha, n)}\\ && +\displaystyle\sum_{i=1}^{c}{\log{\mathcal{N}(\mu_i;\mu_0, (\beta\mathbf{\Lambda}_i)^{-1})\mathcal{W}(\mathbf{\Lambda}_i;\nu, \mathbf{S})}}\\ && +\displaystyle\sum_{i=1}^{c}\sum_{x_k\in\omega_i}{\log{\mathcal{N}(x_k;\mu_i, \mathbf{\Lambda}_i^{-1})}} \end{eqnarray}$

となります。 $\log{\alpha}$ と $\log{\text{AF}(\alpha, n)}$ は4.で事前に計算済みですのでその値を利用できます。この対数事後確率の初期値を計算したらログ用に対数事後確率の最大値 $v_{max}$ と最大値を取った試行回のインデックス $v_{maxi}$ を初期化します。つまり

$\begin{eqnarray} v_{max} &=& \log{v}\\ v_{maxi} &=& 0 \end{eqnarray}$

と初期化します。

学習

概要

学習を行う関数では総学習回数 $iter$ と打ち切り回数 $limit$ を引数として取ることとします。打ち切り回数 $limit$ で指定した回数だけ更新が滞れば収束したものとみなし学習を終了します。また打ち切り回数だけ滞ることがなくとも総学習回数 $iter$ で指定した回数だけ学習を繰り返した場合も学習を終了することとします。ほかに現在の学習回数を表す $iteration$ と現在の停滞回数を表す $stagnation$ の二つを変数として用いることとします。具体的には

$\begin{eqnarray} &&iteration=0\\ &&stagnation=0\\ \\ &&\mathbf{while}\quad stagnation\lt limit\quad\mathbf{and}\quad iteration\lt iter\\ &&\quad\quad\boldsymbol{s}\text{の更新}\\ &&\quad\quad\boldsymbol{\theta}\text{の更新}\\ \\ &&\quad\quad v_{new}=\text{更新後の対数事後確率}\log{p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})}\\ \\ &&\quad\quad \mathbf{if}\quad v_{new}\gt v_{max}\\ &&\quad\quad\quad v_{max} = v_{new}\\ &&\quad\quad\quad v_{maxi} = iteration+1\\ &&\quad\quad\quad stagnation=0\\ &&\quad\quad \mathbf{else}\\ &&\quad\quad\quad stagnation=stagnation+1\\ \\ &&\quad\quad iteration=iteration+1\\ \end{eqnarray}$

のように学習を行います。テキストでは毎回の学習の後に $v_{new}$ が $v_{max}$ を上回っていた場合のみ $\boldsymbol{s}$ に更新を反映させていますが、このやり方だと学習のたびにギブスサンプリングのバーンイン期間の結果を捨てる必要があり計算量が大きくなります。また、もしバーンイン期間を考慮せずにテキスト通りの実装をしてしまうと正しく $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ に従った生成が出来ません。何を初期値に用いるとしても初期値は分布 $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ において確率の低い値でしょう。この状態からただギブスサンプリングを行ってもバーンイン期間がないために生成される $(\boldsymbol{s}, \boldsymbol{\theta})$ は初期値に依存した確率の低い（対数事後確率が小さい）ものになります。こうなるとなかなか $(\boldsymbol{s}, \boldsymbol{\theta})$ が更新されず、更新が起こらないために次のギブスサンプリングでもやはり初期値に依存した確率の低い値が生成され続けるという羽目になります。それを避けるためにバーンイン期間を設けない代わりに $(\boldsymbol{s}, \boldsymbol{\theta})$ を無条件で更新し続け、対数事後確率が最も大きかった試行回を解として出力する形に変更しています。このやり方であれば試行回が進むほど $(\boldsymbol{s}, \boldsymbol{\theta})$ は初期値ではなく $p(\boldsymbol{s}, \boldsymbol{\theta}\vert\boldsymbol{x})$ に従って生成されるようになります。

1. $\boldsymbol{s}$ の更新

$\boldsymbol{s}$ の更新では各 $s_k$ をギブスサンプリングで順次更新していき一通りすべての $s_k$ を更新し終えたところで $\boldsymbol{s}$ の更新を一度やり終えたものとします。具体的には各 $k$ について以下のような処理を繰り返します。

$x_k$ が所属しているクラスタを $\omega_i$ としたとき、そのクラスタから $x_k$ を除きます。具体的には $n_i$ を $-1$ します。
もし $x_k$ を除くことで $\omega_i$ に所属するパターンがなくなる場合、つまり $n_i=0$ となる場合には以下の処理を行います。

・クラスタ数 $c$ を $-1$ する

・ $\boldsymbol{s}_{-k}$ のうちインデックスが $i$ よりも大きいものをすべて $-1$ する

・ $n_i, \mu, \mathbf{\Lambda}$ のうち $i$ 番目のものを削除する

$x_k$ に新しいクラスタを割り当てます。つまり以下を計算して割り当てます。

$\begin{eqnarray}\begin{cases}\displaystyle\frac{n_i'}{\alpha +n-1}\cdot \mathcal{N}(x_k;\mu_i, \mathbf{\Lambda}_i) &(s_k=\omega_i)\\\displaystyle\frac{\alpha}{\alpha +n-1}\cdot\left(\displaystyle\frac{\beta}{(1+\beta)\pi}\right)^{\frac{d}{2}}\cdot\displaystyle\frac{\vert \mathbf{S}_b\vert^{\frac{\nu+1}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1}{2}\right)}{\vert\mathbf{S}\vert^{\frac{\nu}{2}}\cdot\Gamma\left(\displaystyle\frac{\nu+1-d}{2}\right)} &(s_k=\omega_{\text{new}})\end{cases}\end{eqnarray}$

ここで上下ともに第一項の分母 $\alpha +n-1$ は最終的に正規化をすることを考えれば無視できます。上の式について各 $i$ について計算した結果をそれぞれ $r_i$ 、下の式について計算した結果を $r_{new}$ とあらわすならば $x_k$ は

$r_1:r_2:\cdots :r_c:r_{new}$

の比でどのクラスタに割り当てられるかが決まります。

既存クラスタ $\omega_i$ に割り当てられた場合 $n_i$ を $+1$ します。また新規クラスタ $\omega_{new}$ に割り当てられた場合は

・クラスタ数 $c$ を $+1$ する

・ $n_i$ に新規クラスタ分を追加し所属パターン数を $1$ とする

・ $x_k$ のみを要素とするクラスタとして $\theta_{new}$ を生成する

・ $\mu, \mathbf{\Lambda}$ に新規クラスタ分を追加し、生成した $\theta_{new}$ を割り当てる

と処理します。 $\theta_{new}$ の生成方法は後述の2. $\boldsymbol{\theta}$ の更新と同じです。

2. $\boldsymbol{\theta}$ の更新

各クラスタ毎に

$\begin{eqnarray} \mathbf{\Lambda}_i &\sim& \mathcal{W}(\nu_c, \mathbf{S}_q)\\ \mu_i &\sim& \mathcal{N}(\mu_c, \mathbf{\Lambda}_c^{-1}) \end{eqnarray}$

のように $\mathbf{\Lambda}_i$ と $\mu_i$ を確率的に生成し更新します。ここで $\nu_c, \mathbf{S}_q, \mu_c, \mathbf{\Lambda}_c$ は

$\begin{eqnarray} \bar{x} &=& \displaystyle\frac{1}{n_i}\displaystyle\sum_{x_k\in\omega_i}{x_k}\\ \nu_c &=& \nu+n_i\\ \mathbf{S}_q^{-1} &=& \mathbf{S}^{-1}+\displaystyle\sum_{x_k\in\omega_i}{(x_k-\bar{x})(x_k-\bar{x})^t}+\displaystyle\frac{n_i\beta}{n_i+\beta}(\bar{x}-\mu_0)(\bar{x}-\mu_0)^t\\ \mu_c &=& \displaystyle\frac{n_i\bar{x}+\beta\mu_0}{n_i+\beta}\\ \boldsymbol{\Lambda}_c &=& (n_i+\beta)\boldsymbol{\Lambda}_i\\ \end{eqnarray}$

となります。

3. 対数事後確率の計算

計算の内容は初期化の「5. 対数事後確率の初期値の計算」と同じです。つまり

$\begin{eqnarray} \log{v} &=& c\log{\alpha}+\displaystyle\sum_{i=1}^{c}{\sum_{l=1}^{n_i-1}{\log{l}}}-\log{\text{AF}(\alpha, n)}\\ && +\displaystyle\sum_{i=1}^{c}{\log{\mathcal{N}(\mu_i;\mu_0, (\beta\mathbf{\Lambda}_i)^{-1})\mathcal{W}(\mathbf{\Lambda}_i;\nu, \mathbf{S})}}\\ && +\displaystyle\sum_{i=1}^{c}\sum_{x_k\in\omega_i}{\log{\mathcal{N}(x_k;\mu_i, \mathbf{\Lambda}_i^{-1})}} \end{eqnarray}$