2014年11月24日月曜日

1次元ロジスティック回帰と確率的勾配降下法


はじめに


 ここでは、1次元ロジスティック回帰と、これを解く際に用いられる確率的勾配降下法についてまとめる。

1次元ロジスティック回帰


 2分類問題を考える。$N$ 個の観測値 $\vec{x}^{\;i},i=1,\cdots,N$ が与えられ、それぞれ0か1のラベルを持つとする。このとき、観測値が従う確率分布モデル $p(c=0,1\;|\;\vec{x}\;;\;\vec{\omega},b)$ を求めたい。ここで、$\vec{\omega},b$ はモデルを調節するパラメータである。各観測値が他の観測値の影響を受けないと仮定すれば、全ての観測値の確率は 次式で与えられる。 \begin{equation} L(\vec{\omega},b)=\prod_{i=1}^{N}\;p(c^{i}\;|\;\vec{x}^{\;i}\;;\;\vec{\omega},b) \label{likelihood} \end{equation} これが最大となるように $(\vec{\omega},b)$ を調節すればよい。計算を容易にするため次式の最小化を考える。 \begin{equation} L(\vec{\omega},b)=-\sum_{i=1}^{N}\;\ln{p(c^{i}\;|\;\vec{x}^{\;i}\;;\;\vec{\omega},b)} \label{nll} \end{equation} 上式を負の対数尤度(Negative Log Likelihood:NLL)関数と呼ぶ。ところで、観測値のラベル $c$ を0か1に制限したので、確率 $p$ は以下のように書くことができる。 \begin{equation} p(c\;|\;\vec{x}\;;\;\vec{\omega},b)=s(1\;|\;\vec{x}\;;\;\vec{\omega},b)^c\;\left\{1-s(1\;|\;\vec{x}\;;\;\vec{\omega},b)\right\}^{1-c} \end{equation} ここで、$s$ は $c=1$ となる確率である。$s$ と $\vec{x}$ を結びつけるモデルが必要である。1次元ロジスティック回帰では次のモデルを考える。 \begin{eqnarray} s(1\;|\;\vec{x}\;;\;\vec{\omega},b) &=& \sigma\left(f(\vec{x}\;;\;\vec{\omega},b)\right)\nonumber \\ \sigma(a) &=& \frac{1}{1+\exp{(-a)}}\label{sigmoid}\\ f(\vec{x}\;;\;\vec{\omega},b) &=& \vec{\omega}^{\;T}\cdot\vec{x}+b \label{projection} \end{eqnarray} 関数 $\sigma$ はシグモイド関数と呼ばれ、以下のような形をしている。
ベクトル $\vec{x}$ を式(\ref{projection})により1次元に射影し、その結果を式(\ref{sigmoid})を使って区間(0,1)に射影し、これを確率とみなすわけである。式(\ref{nll})をパラメータ $(\vec{\omega},b)$ で微分するのであるが、その前に少しだけ式を一般化する。いま、$\vec{x}$ に1次元加えたベクトルを導入する。 \begin{equation} \vec{X}=(\vec{x},1) \end{equation} 最後の次元の値は常に1である。さらに、$\vec{\Omega}$ を次式で定義する。 \begin{equation} \vec{\Omega}=(\vec{\omega},b) \end{equation} このとき関数 $f$ は \begin{equation} f(\vec{X}\;;\;\vec{\Omega}) = \vec{\Omega}^{\;T}\cdot\vec{X} \label{new_projection} \end{equation} と書くことができる。この表式を使うと $L$ は \begin{equation} L(\vec{\Omega})=-\sum_{i=1}^{N} \left\{ c^{i}\ln{s(1\;|\;\vec{X}^{\;i}\;;\;\vec{\Omega})} + (1-c^{i}) \ln{ \left( 1-s(1\;|\;\vec{X}^{\;i}\;;\;\vec{\Omega}) \right) } \right\} \end{equation} となる。ただし、 \begin{eqnarray} s(1\;|\;\vec{X}\;;\;\vec{\Omega}) &=& \sigma\left(f(\vec{X}\;;\;\vec{\Omega})\right)\nonumber \\ \sigma(a) &=& \frac{1}{1+\exp{(-a)}}\label{new_sigmoid}\\ f(\vec{X}\;;\;\vec{\Omega}) &=& \vec{\Omega}^{\;T}\cdot\vec{X} \label{new_new_projection} \end{eqnarray} である。この式を $\vec{\Omega}$ で微分する。シグモイド関数なので計算は容易である。 \begin{equation} \frac{\partial L(\vec{\Omega})}{\partial \vec{\Omega}}=-\sum_{i=1}^{N}\;\vec{X}^{\;i} \left(c^{i}-s(1\;|\;\vec{X}^{\;i}\;;\;\vec{\Omega})\right) \label{gradient} \end{equation} これを0と置いたときの解が求める $\vec{\Omega}$ であるが、これ以上の解析計算は不可能である。そこで、勾配降下法を使って数値計算を行なう。手順は以下の通りである。
  1. $\vec{\Omega}$に適当な初期値を与える。
  2. 次式で $\vec{\Omega}$ を更新する。 \begin{equation} \vec{\Omega}\leftarrow \vec{\Omega} - \eta \frac{\partial L(\vec{\Omega})}{\partial \vec{\Omega}} \end{equation} ここで $\eta$ は学習率と呼ばれる正の量である。収束速度を調節するパラメータである。
  3. $L$ の変化量がしきい値以下となるまで上記を繰り返す。あるいは、エラー $E$ \begin{equation} E=\sum_{i=1}^{N}\;\left|c^{i}-s(1\;|\;\vec{X}^{\;i}\;;\;\vec{\Omega})\right| \end{equation} の変化量がしきい値以下となるまで繰り返しても良い。この量は理想的には0になって欲しい量である。

確率的勾配降下法(Stochastic Gradient Descent:SGD)


 式(\ref{gradient})を見て分る通り、通常の勾配降下法では $\vec{\Omega}$ の更新ごとに全観測値についての和を取る必要がある。大規模データの場合、これでは大変時間がかかる。そこで、$\vec{\Omega}$ の更新に使うデータをランダムに取り出すことを考える。この手法を確率的勾配降下法と呼ぶ。 手順は以下の通りである。
  1. データをシャッフルしたあと、複数のグループに分割する。たとえば、1000個のデータを10個のデータのグループに分割する。各グループはミニバッチ(minibatch)と呼ばれる。
  2. $\vec{\Omega}$に適当な初期値を与える。
  3. 次式で $\vec{\Omega}$ を更新する。
    • $\mbox{for}\;\;\;\;k = 1, \cdots,$ \begin{equation} \vec{\Omega}\leftarrow \vec{\Omega} - \eta \left\{ -\frac{1}{|\Lambda^{k}|} \sum_{i\;\in\;\Lambda^k}\;\vec{X}^{\;i} \left(c^{i}-s(1\;|\;\vec{X}^{\;i}\;;\;\vec{\Omega})\right) \right\} \label{grad} \end{equation}
    ここで $\Lambda^{k}$ は $k$ 番目のミニバッチを、$|\Lambda^{k}|$ はそこに含まれるデータの数を表す。
  4. 先の勾配法と同じように何らかの量がしきい値以下となるまで上記を繰り返す。

実装


 2次元平面に、(5,5)と(0,0)を中心にした2つの正規分布を使って1000個ずつ点を生成し、2分類問題を行なった。ミニバッチのサイズを10、学習率を0.1とした。Theanoを使って実装を行なった。
  • 45,46行目:行列とベクトルのシンボルを作成する。前者は座標を後者はラベルを格納する。
  • 49行目:パラメータ $\vec{\Omega}$ を表す共有変数である。
  • 55行目:$f$ を表すTheano Expression Graphである。
  • 58行目:$\sigma$ を表すTheano Expression Graphである。
  • 61行目:Negative Log Likelihood(NLL) を表すTheano Expression Graphである。
  • 64行目:ミニバッチ内でのNLLの平均値である。
  • 67行目:$\vec{\Omega}$についての微分を表すTheano Expression Graphである。
  • 70行目:学習率 $\eta$ を表す共有変数である。
  • 71行目:$\vec{\Omega}$ の更新を表現するTheano Functionである。
  • 76行目:$\eta$ の更新を表現するTheano Functionである。
  • 79行目:サンプルを作成する。
  • 84〜92行目:勾配降下法を実行する
  • 下図における青線が、識別関数 $f(\vec{X}\;;\;\vec{\Omega}) = \vec{\Omega}^{\;T}\cdot\vec{X}=0$ である。
    このときのエラー $E$(プログラム中のリストerrors)を下に示す。合計2000個のデータを200個のミニバッチに分割し、これを10回繰り返したので横軸の最大値は2000である。 振動しながら小さくなっていく様子が分かる。

    2014年11月3日月曜日

    Restricted Boltzmann Machine


    はじめに


     以前、Boltzmann Machineの基礎を解説し、その最後に隠れ変数の意味について述べた。今回はRestricted Boltzmann Machineについて解説する。

    Restricted Boltzmann Machine


     Restricted Boltzmann Machineとは、可視層(入力層)と隠れ層(出力層)の2層から構成されるBoltzmann Machineである。ユニット間の結合は層間にのみ存在し、同一層内には存在しない(下図参照)。これが名前の由来である。
    図の $\Gamma$ は可視ユニットの集合を、$\Lambda$ は隠れユニットの集合を表す。 エネルギー関数は次式で定義される。 \begin{equation} \Phi({\bf V},{\bf H}|{\boldsymbol \Theta}) = -\sum_{i\;\in\;\Gamma}\; \theta_{i}^V\;V_{i} -\sum_{j\;\in\;\Lambda}\;\theta_{j}^H\;H_{j} - \sum_{i\;\in\;\Gamma}\;\sum_{j\;\in\;\Lambda} \omega_{ij}\;V_{i}\;H_{j} \label{energy} \end{equation} ただし、${\boldsymbol \Theta}=({\boldsymbol \theta}^V, {\boldsymbol \theta}^H, {\boldsymbol \omega})$ とした。 可視層と隠れ層のバイアス項は区別される。 ${\bf V}$ と ${\bf H}$ はそれぞれ可視確率変数、隠れ確率変数を表す。

    学習方程式


     Restricted Boltzmann Machineの学習モデルは先の解説と同様にBoltzmann関数である。 \begin{equation} P_{B}({\bf V}, {\bf H}|{\boldsymbol \Theta}) = \frac{1}{ Z_{B}({\boldsymbol \Theta}) } \exp\left[-\Phi({\bf V},{\bf H}|{\boldsymbol \Theta})\right] \label{boltzmann-1} \end{equation} 式を扱い易くするため、統計物理学で使用される自由エネルギー $F({\boldsymbol V}|{\boldsymbol \Theta})$ \begin{equation} F({\boldsymbol V}|{\boldsymbol \Theta})= -\ln{ \left( \sum_{H_1=\pm1}\cdots\sum_{H_m\pm1}\; \exp \left[ -\Phi({\bf V},{\bf H}|{\boldsymbol \Theta}) \right] \right) } \label{free} \end{equation} を導入する。このとき正規化因子(分配関数)は次式で与えられる。 \begin{equation} Z_{B}({\boldsymbol \Theta})= \sum_{V_1=\pm1}\cdots\sum_{V_n=\pm1}\; \exp \left[ -F({\boldsymbol V}|{\boldsymbol \Theta}) \right] \end{equation} 式(\ref{energy})を式(\ref{free})に代入する。 \begin{eqnarray} F({\boldsymbol V}|{\boldsymbol \Theta}) &=& -\ln{\left( \sum_{H_1}\cdots\sum_{H_m}\; \exp \left[ \sum_{i\;\in\;\Gamma}\; \theta_{i}^V\;V_{i} +\sum_{j\;\in\;\Lambda}\;\theta_{j}^H\;H_{j} +\sum_{i\;\in\;\Gamma}\;\sum_{j\;\in\;\Lambda} \omega_{ij}\;V_{i}\;H_{j} \right] \right) } \nonumber \\ &=& -\ln{\left( \exp \left[ \sum_{i\;\in\;\Gamma}\; \theta_{i}^V\;V_{i} \right] \sum_{H_1=\pm 1}\cdots\sum_{H_m=\pm 1} \;\exp \left[ \sum_{j\;\in\;\Lambda}\;\theta_{j}^H\;H_{j} +\sum_{i\;\in\;\Gamma}\;\sum_{j\;\in\;\Lambda} \omega_{ij}\;V_{i}\;H_{j} \right] \right)} \nonumber \\ &=& -\ln{\left( \exp \left[ \sum_{i\;\in\;\Gamma}\; \theta_{i}^V\;V_{i} \right] \sum_{H_1=\pm 1} \;\exp \left[ \left\{\theta_{1}^H+\sum_{i\;\in;\Gamma}\;\omega_{i1}\;V_{i}\right\}H_1 \right] \cdots \sum_{H_m=\pm 1} \;\exp \left[ \left\{\theta_{m}^H+\sum_{i;\in\;\Gamma}\;\omega_{im}\;V_{i}\right\}H_m \right] \right)} \nonumber \\ &=& -\ln{\left( \exp \left[ \sum_{i\;\in\;\Gamma}\; \theta_{i}^V\;V_{i} \right] \prod_{j\;\in\;\Lambda}\;2\;\cosh{\left\{\theta_{j}^H+\sum_{i\;\in\;\Gamma}\;\omega_{ij}\;V_{i}\right\}} \right)} \nonumber \\ &=& -\sum_{i\;\in\;\Gamma}\;\theta_i^V\;V_i - \sum_{j\;\in\;\Lambda}\;\ln{ \left( 2\cosh{\left\{\theta_j^H+\sum_{i\;\in\;\Gamma }\omega_{ij}\;V_i\right\}} \right) } \label{free2} \end{eqnarray} $P_{B}({\bf V}, {\bf H}|{\boldsymbol \Theta})$を隠れ変数について周辺化したものは、自由エネルギー $F$ を使って 書くことができる。 \begin{eqnarray} P_{V}({\bf V}|{\boldsymbol \Theta}) &\equiv& \sum_{H_1}\cdots\sum_{H_m}\;P_{B}({\bf V}, {\bf H}|{\boldsymbol \Theta}) \nonumber\\ &=&\frac{1}{Z_B({\boldsymbol \Theta})}\; \exp{ \left[ -F\left({\boldsymbol V}|{\boldsymbol \Theta}\right) \right] } \label{marginalv} \end{eqnarray} この式を使って次式を最大化することを考える。 \begin{equation} L_{D}({\boldsymbol \Theta})=\sum_{\mu=1}^{N}\;\ln{P_{V}(\vec{v}^{\;\mu}|{\boldsymbol \Theta})} \label{logld} \end{equation} ここで、$\vec{v}^{\;\mu}$は観測値である。 議論を一般化するため、パラメータ群 ${\boldsymbol \Theta}$ を代表するパラメータ $\xi$ を考え、これによる偏微分を 考える。 \begin{equation} \frac{\partial L_D}{\partial \xi} = 0 \Leftrightarrow E_V\left(\frac{\partial F({\boldsymbol V}|{\boldsymbol \Theta})}{\partial \xi} \right) = \frac{1}{N}\;\sum_{\mu}\; \frac{\partial F(\vec{v}^{\;\mu}|{\boldsymbol \Theta})}{\partial \xi} \label{general} \end{equation} ここで、${\bf V}$ についての期待値 \begin{equation} E_V(\nu)=\sum_{V_1}\cdots \sum_{V_n}\; \nu\;P_{V}({\bf V}|{\boldsymbol \Theta}) \end{equation} を定義した。 式(\ref{general})の左辺は確率変数についての期待値、右辺は観測値ついての期待値である。 $\xi=\{\theta_i^V, \theta_j^H, \omega_{ij}\}$ として次式を得る。 \begin{eqnarray} \frac{\partial L_D}{\partial \theta_i^V} = 0 &\Leftrightarrow& E_V(V_i) = \frac{1}{N}\;\sum_{\mu}\;v_i^{\;\mu} \label{le-1}\\ \frac{\partial L_D}{\partial \theta_j^H} = 0 &\Leftrightarrow& E_V \left( \tanh{ \lambda_j^{H} } \right) = \frac{1}{N}\;\sum_{\mu}\; \tanh{\lambda_j^{H,\;\mu}} \label{le-2}\\ \frac{\partial L_D}{\partial \omega_{ij}} = 0 &\Leftrightarrow& E_V\left( V_i \tanh{ \lambda_j^{H} } \right) = \frac{1}{N}\;\sum_{\mu}\; v_i^{\;\mu} \tanh{\lambda_j^{H,\;\mu}} \label{le-3} \end{eqnarray} ここで \begin{eqnarray} \lambda_j^{H,\;\mu}&=&\theta_j^H+\sum_{k\;\in\;\Gamma}\omega_{kj}\;v_k^{\;\mu} \\ \lambda_j^{H}&=&\theta_j^H+\sum_{k\;\in\;\Gamma}\omega_{kj}\;V_k \end{eqnarray} とした。以下に $y=\tanh(x)$ のグラフを示す。
    式(\ref{le-1})(\ref{le-2})(\ref{le-3})が Restricted Boltzmann Machine の学習方程式である。3式とも、左辺は確率変数 ${\bf V}$ についての期待値であり、「計算量の爆発」問題を抱えている。近似的解法が必要である。

    RBMの性質


     式(\ref{marginalv})とベイスの公式を使うと、確率変数 ${\bf V}$ を固定したもとで隠れ変数 ${\bf H}$ が$\pm 1$になる確率を求めることができる。 \begin{eqnarray} P_{H|V}({\bf H}=\pm 1|{\bf V},\Theta)&=&\frac{P_B({\bf V},{\bf H}|\Theta)}{P_V({\bf V}|\Theta)} \nonumber \\ &=&\prod_{j\;\in\;\Lambda}\;\frac{1}{1+\exp{\left(\mp\;2\;\lambda_{j}^{H}\right)}} \label{prior-1} \end{eqnarray} また、$P_B({\bf V},{\bf H}|\Theta)$ に対し、可視変数で周辺化し同様な計算を繰り返すと、確率変数 ${\bf H}$を固定したもとで可視変数 ${\bf V}$ が$\pm 1$になる確率を求めることができる。 \begin{equation} P_{V|H}({\bf V}=\pm 1|{\bf H},\Theta)=\prod_{i\;\in\;\Gamma}\;\frac{1}{1+\exp{\left(\mp\;2\;\lambda_{i}^{V}\right)}} \label{prior-2} \end{equation} ここで \begin{eqnarray} \lambda_j^H&=&\theta_j^H+\sum_{k\;\in\;\Gamma }\omega_{kj}\;V_k \\ \lambda_i^V&=&\theta_i^V+\sum_{k\;\in\;\Lambda}\omega_{ik}\;H_k \end{eqnarray} とした。以下に $y=1/(1+\exp{(-x)})$ のグラフを示す。
    0から1の間を変化する関数、すなわち確率であることが確認できる。 式(\ref{prior-1})(\ref{prior-2})はいずれも、ユニットごとの確率の積となっている。これは、Restricted Boltzmann Machine に於いて、層内結合が存在しないためである。

     ところで、先の解説で、隠れ変数の導入はモデルの表現能力を高めると述べた。 式(\ref{marginalv})の形から、自由エネルギー $F$ は$P_B$における$\Phi$ に相当するものであることが分る。 $F$の第二項に対し、$\omega_{ij}$ を微小量としてテイラー展開を行なう。 \begin{eqnarray} \ln{2\cosh{\left(\theta_j^H+\sum_{k\;\in\;\Gamma }\omega_{kj}\;V_k\right)}}&=& \ln{2\cosh{\theta_j^H}}+\tanh{\theta_j^H} \cdot \sum_{k\;\in\;\Gamma}\omega_{kj}\;V_k \\ &+&\frac{1}{2} {\rm sech}^2 \theta_j^H \cdot \left( \sum_{k\;\in\;\Gamma}\omega_{kj}\;V_k \right)^2 \\ &-&\frac{1}{3} {\rm sech}^2 \theta_j^H \cdot \tanh{\theta_j^H} \cdot \left( \sum_{k\;\in\;\Gamma}\omega_{kj}\;V_k \right)^3\\ &+&{\it O}(\omega^4) \end{eqnarray} 最初に導入したエネルギー関数(\ref{energy})は可視ユニット間の相互作用を含まないが、"有効"エネルギー関数にはこの相互作用が現れている。これは、可視ユニット同士が隠れユニットを介して結合することを意味する。相互作用の出現によりモデルの表現能力が高まるのである。

    近似的解法(Contrastive Divergence)


     ここでは、確率変数についての期待値 $E_V(\nu)$ \begin{equation} E_V(\nu)=\sum_{V_1=\pm1}\cdots \sum_{V_n=\pm1}\; \nu(V_1,\cdots,V_n|{\boldsymbol \Theta})\;P_{V}(V_1,\cdots,V_n|{\boldsymbol \Theta}) \end{equation} を評価する近似法について述べる。右辺は $2^{n}$ 個の項の和であるが、これを観測値と同じ個数($N$)の項の和で近似する。 \begin{equation} E_V(\nu)\sim \frac{1}{N}\sum_{\mu=1}^{N}\; \nu(V_1^{\mu},\cdots,V_n^{\mu}|{\boldsymbol \Theta}) \label{app} \end{equation} 話しを進める前に、式(\ref{prior-1})(\ref{prior-2})から各ユニットの条件付き確率を求めておく。 \begin{eqnarray} P_{H|V} (H_j=\pm 1|{\bf V},\Theta) &=&\frac{1}{1+\exp{\left(\mp\;2\;\lambda_{j}^{H}\right)}},\;\; j\;\in\;\Lambda \label{unit1} \\ P_{V|H} (V_i=\pm 1|{\bf H},\Theta) &=&\frac{1}{1+\exp{\left(\mp\;2\;\lambda_{i}^{V}\right)}},\;\; i\;\in\;\Gamma \label{unit2} \end{eqnarray} まず最初に、${\bf V}^{\mu}=(V_1^{\mu},\cdots,V_n^{\mu})$ の初期値 ${\bf V}_0^{\;\mu}$ として観測値 $\vec{v}^{\;\mu}$ を採用する。 さらに、確率(\ref{unit1})を使って、${\bf H}^{\mu}=(H_1^{\mu},\cdots,H_m^{\mu})$ の初期値 ${\bf H}_0^{\;\mu}$ を計算する。 ${\bf H}^{\mu}_0$ と確率(\ref{unit2})を使って ${\bf V}^{\mu}_1$ を計算する。これを $T$ 回繰り返す。 \begin{equation} {\bf V}^{\mu}_0 \rightarrow {\bf H}^{\mu}_0 \rightarrow {\bf V}^{\mu}_1 \rightarrow {\bf H}^{\mu}_1 \cdots {\bf H}^{\mu}_{T-1} \rightarrow {\bf V}^{\mu}_T \nonumber \end{equation} ${\bf V}^{\mu}_T$ を式(\ref{app}) に代入する。 \begin{equation} E_V(\nu)\sim \frac{1}{N} \sum_{\mu=1}^{N}\; \nu({\bf V}^{\mu}_T|{\boldsymbol \Theta}) \end{equation} $T=1$ でも精度よく計算できることが知られている。

    参考文献

    1. ディープボルツマンマシン入門 ボルツマンマシン学習の基礎, 安田宗樹, 人工知能学会誌 28巻 3号, 2013年5月.
    2. Restricted Boltzmann Machines (RBM)
    3. コンピュータビジョン最先端ガイド6, 第4章 ディープラーニング