nykergoto’s blog

機械学習とpythonをメインに、雑多な内容をとりとめなく扱うブログです。

ニューラルネットへのベイズ推定 - Bayesian Neural Network

ニューラルネットワーク過学習防止としてもちいられる Dropout がネットワークの重みに対してベイズ推定を行っていることに相当する、ということを述べた論文について紹介します。

くわえてそれを実装して簡単な実験をやってみます。

Ref

ベイズ推定

はじめにベイズ推定について簡単に。
ベイズ推定で目指す目的とは何でしょうか。

仮定として, $N$ 個の入力データ $X=\{x_1, \cdots. x_N \}$ と、それに対応する出力 $Y=\{ y_1, \cdots, y_N \}$ をすでに観測しているとします。またデータは独立同分布から生成されている, とします。
この状態で、新しい入力 $x$ が入ってきたときに出力 $y$ を予測するという問題を考えます。
ベイズ推定で求めたいのは、データが与えられたときの、隠れ変数の事後分布 $p(w \mid X, Y)$ です。

重みの事後分布さえ入手できてしまえば、それを用いて新しいデータ $\hat{x}$ が与えられたとき、その予測値 $\hat{y}$ が従う分布は以下で計算することが可能です。

$$p(\hat{y} \mid \hat{x}, X, Y) = \int p(\hat{y} \mid w, \hat{x}) p(w \mid X, Y) dw$$

ですから、ベイズ推定の立場にたてば、事後分布 $p(w \mid X, Y)$ をいかにして求めていくか, ということが問題となります。

事後分布の計算

直接事後分布が計算できるような単純なモデルでは、先のアイディアに基づいて事後分布を計算すれば良いです。

しかし、複雑なモデルでは、陽に事後分布が計算できない場合が出てきます。 このような仮定のもとでは、何らかの関数で事後分布を近似する必要があります。
以下ではあるパラメータ $\theta$ を持つ関数 $p_\theta(w)$ を事後分布に近づいけていく, という方法を考えます。
このとき、2つの分布の近さを測る距離として, 2つの確率分布から実数空間へ射影する関数として KL-Divergence (以下では ${\rm KL}$ と表記) を用います。
KL距離は距離空間ではありませんが, 確率分布 $p,q$ に対して, ${\rm KL}(q, p) = 0$ のとき $q=p$ が成立します。よって KL 距離の意味で最小化を行えば, 分布 $q_\theta$ は $p(w\mid X, Y)$ に近づく事が期待されます.
最小化する関数を変形すると、以下のようになります。

$$\begin{align}{\rm KL}(q_{\theta}(w) | p(w | X, Y))&=\int q_{\theta}(w) \log \left( \frac{q_{\theta}(w)}{p(w | X,Y)} \right)dw \\&=\int q_{\theta}(w) \log \left( q_{\theta}(w) \frac{p(X,Y)}{p(X,Y|w) p(w)} \right) dw \\&=\int q_{\theta}(w) \left( \log \frac{q_{\theta}(w)}{ p(w)} + \log \frac{p(Y|X) p(X)}{p(Y|X, w) p(X|w)} \right) dw \\&=\int q_{\theta}(w) \left( \log \frac{q_{\theta}(w)}{ p(w)} + \log \frac{p(Y|X)}{p(Y|X, w)} \right) dw \\&= \int q_{\theta}(w) \log p(Y|X) dw - \int q_{\theta}(w) \log p(Y|X,w)dw + {\rm KL} (q_{\theta}(w) | p(w)) \\&= \log p(Y|X) - \int q_{\theta}(w) \log p(Y| X, w)dw + {\rm KL} (q_{\theta}(w) | p(w))\end{align}$$

ここで第4行目の変形で $p(X|w)=p(X)$ を用いました。今 Variance Inference を ${\rm VI}$ と表し

$${\rm VI}(\theta) = - \int q_{\theta}(w) \log p(Y|X,w)dw + {\rm KL} (q_{\theta}(w) | p(w))$$

と定義します。すると上記の式は以下のように変形することができます.

$${\rm KL}(q_{\theta}(w) | p(w | X, Y)) = {\rm VI}(\theta) + \log p(X,Y)$$

ここで右辺最終項 $\log p(X,Y)$ は定数ですから, 結局 KL 距離の最小化問題は以下のように書き換えられます。

$$\min_{\theta} {\rm KL}(q_{\theta}(w) | p(w | X, Y)) = \min_{\theta}{\rm VI}(\theta)$$

よって ${\rm VI}$ を $\theta$ に関して最小化することは, KL-Divergence の意味での最小化と一致する事がわかります。
また ${\rm VI}$ の第一項

$$-\int q_{\theta}(w) \log p(Y|X,w)dw$$

の部分は, $0 \le p(Y|X,w) \le 1$ であることを考慮すると必ず正の値です。また得られたデータ $Y$ が現れる確率 $p(Y|X,w)$ の値が高い $w$ に対して、分布 $q_\theta(w)$ も高い値を取るときに小さくなり、データへの当てはまりを表していることがわかります。

一方第二項の

$${\rm KL} (q_{\theta}(w) | p(w))$$

については, 事前分布 $p(w)$ と事後分布の推定密度関数 $q_\theta(w)$ との距離を表しています。これは事前分布からかけ離れた分布へのペナルティを与える項に相当しており、一般に正則化項と呼ばれるものです。(簡単な例で言えば互いにガウス分布であるときこのKLダイバージェンス部分は $L2$ ノルムとなり線形回帰であれば Ridge 回帰とよばれる枠組みになります)
このように自分が提案する分布 $q$ をもとめたい分布 $p$ に近づけていく方法を変分推論 (Variance Inference) とよびます.

ニューラルネットワークに対する既存の VI

つぎにニューラルネットワークに対して Variance Inference を行う既存の枠組みについて振り返ってみます。近年のニューラルネットワークへのベイズ推定のアプローチの一つである Hinston and Van Camp, 1993 では近似する分布 $q_\theta$ にたいして, すべての重みに対して独立したガウス分布を仮定しました。

すなわち $I$ 層の隠れ層を持ち, その各 $i$ 層で重み $w_{ijk}$ を持つネットワークに対して $m_{ijk}, \sigma_{ijk} \in \mathbb{R}$ とおいたとき

$$q_\theta (w) = \prod_{i, j ,k} N(w_{ijk}\mid m_{ikj}, \sigma_{ijk}^2)$$

と記述することになります。 この分布を仮定した場合の最適化はとてもむずかしく, 論文中では一つの隠れ層を持つニューラルネットワークに対しての実験にとどまっていました。 実際理論的には事後分布の推論にガウス分布を用いると一層の隠れ層をもつニューラルネットに対しては数式上綺麗に解析が行えるものの、実験上性能がよくありませんでした。理由としては、この方法が本来重要である重み同士の関係性を記述できないことなどが挙げられています。

これに対し Barver and Bishop, 1998 では重みに対して混合ガウスを仮定したモデルを提案しました。これによって同一層の重みに対して関連性を考慮することが可能となりました。しかし、その反面計算量が増大してしまうため、複雑なモデル等に対してうまく働きません。

Variance Inference の近似

そのままの形ではなかなか良い結果を得ることができていなかった変分推論ですが、今回は期待値の意味で一致する近似式を用いて最適化することを考えていきます。


まず $X,Y$ は独立であると仮定しているので $p(Y | w,X) = \prod_{i} p(y_i | w, x_i )$ が成り立ちます。これを用いて Variance Inference を変形すると以下のようになります。

$$\begin{align*}{\rm VI}(\theta) &= - \int q_{\theta}(w) \log p(X,Y|w)dw + {\rm KL} (q_{\theta}(w) | p(w)) \\&= -\sum_{i=1}^N \int q_{\theta}(w) \log p(y_i | x_i, w) dw + {\rm KL} (q_{\theta}(w) | p(w)) \\&= -\sum_{i=1}^N \int q_{\theta}(w) \log p(y_i | f^w(x_i)) dw + {\rm KL} (q_{\theta}(w) | p(w))\end{align*}$$

ここで $f^w(x_i)$ は重み $w$ を持つモデルに $x_i$ が入力されたときの出力を表します。

この形式を取り扱う上での難しさは, 主に以下の二点にあります

  1. $w$ に関する積分の部分が扱いやすい形式ではない
  2. データ $N$ に対する和を取らなくてはならないため, データの数が多くなると計算が困難になる

この内 2 に関してはすべてのデータのうちで \( M \) 個だけサンプルする (mini-batchを計算する) ことで回避することができます. この近似を行った $\hat{{\rm VI}}$ は以下のようになります。

$$\hat{{\rm VI}} = - \frac{N}{M} \sum_{i \in S} \int q_{\theta}(w) \log p(y_i | f^w(x_i)) dw + {\rm KL} (q_{\theta}(w) | p(w))$$

ここで, $S$ はサンプルされた添字の集合を表します。

しかし理由 1 によりこの計算は難しいままです。

$w$ の積分の近似

$w$ に関する積分を近似することを考えます。 積分の部分を見ると, $q_\theta$ と対数尤度との掛け算になっています。よって $q_\theta$ から $w$ をサンプリングすることができれば積分計算の近似を行うことができます。

しかし, $q_\theta$ の形式には仮定が置かれておらず任意の分布を取ることが可能です。したがってこの分布からサンプリングすることはできません。

そこで $q_\theta(w)$ がパラメータを持たない別の分布 $p(\epsilon)$ を用いて, $w = g(\theta, \epsilon)$ で表現できるという仮定を加えてみます。そうすると分布 $q$ に対する情報がわかっていなくても、分布 $p(\epsilon)$ からサンプルした値 $\hat{\epsilon}$ を用いて提案分布からのサンプル $\hat{w} = g(\theta, \hat{\epsilon})$ を生成することが可能となります。

平均 0 分散 1 のガウス分布 (正規分布) のサンプル値から 任意の分散と平均値をもつガウス分布のサンプルを作成することなどがこれに相当します。

これを用いて先の VI の式中の $q_\theta$ を $p(\epsilon)$ で表現することで, 分布の部分から $\theta$ を取り除きます。

$$ \begin{align} \hat{{\rm VI}} &= - \frac{N}{M} \sum_{i \in S} \int q_{\theta}(w) \log p(y_i | f^w(x_i)) dw + {\rm KL} (q_{\theta}(w) | p(w)) \\ &= - \frac{N}{M} \sum_{i \in S} \int p(\epsilon) \log p(y_i | f^{g(\theta, \epsilon)}(x_i)) d\epsilon + {\rm KL} (q_{\theta}(w) | p(w)) \end{align} $$

こうなると積分はパラメータを持たない分布 $p(\epsilon)$ からのサンプリングを行うことで効率的に近似をおこなうことが可能になります。 サンプリングをミニバッチのサンプルと同時に行うとし、サンプルされた実現値を $\epsilon_i \sim p(\epsilon) (i = 1, \cdots, M)$ とします。
するとこのモンテカルロ法によって積分が近似された $\hat{{\rm VI}}_{MC}$ は以下のようになります。

$$\hat{{\rm VI}}_{MC} = - \frac{N}{M} \sum_{i \in S} \log p(y_i | f^{g(\theta, \epsilon_i)}(x_i)) + {\rm KL} (q_{\theta}(w) | p(w))$$

このとき $\epsilon, S$ に対する期待値を取ると, $\mathbb{E}_{S, \epsilon} \left[\hat{{\rm VI}}_{MC} \right] = {\rm VI}$ が成立します。

この新しい $\hat{{\rm VI}}_{MC}$ を目的関数として例えば勾配法を用いて $\theta$ について最適化を行えば、期待値を取れば元の Variance Inference を最適化することと同値です。よって最適化の各 $t$ ステップにおいて

$$ \displaystyle \theta_{t} = \theta_{t-1} + \eta \frac{\partial}{\partial\theta} \hat{{\rm VI}}_{MC} $$

によりパラメータ $\theta$ を更新すれば良いことがわかります。

Dropout による学習

ここで一旦 Variance Inference のことをおいておいて, ニューラルネットワーク過学習を抑える手法の一つである Dorpout について考えます。
Dropout とは訓練データが与えられたとき, 各層においてすべての隠れノードを用いて出力を行わず, ランダムに選ばれたノードの値のみを用いて出力をし, backword においても出力に関わったノードの値のみを更新する, という方法です。

単純に一層の隠れ層を持つネットワークを考え, 入力から隠れ層, 隠れ層から出力層への重みをそれぞれ $M_1 \in \mathbb{R}^{n\times m}, M_2 \in \mathbb{R}^{m\times l} $とします.
また隠れ層の定数 $b \in \mathbb{R}^{m}$, 活性化関数 $\sigma$ とします。

今入力として $x \in \mathbb{R}^n$ のデータが与えられたとします。この時 $n$ 次元上の確率 $p_1 (0\le p_1\le 1)$ のベルヌーイ分布 $Q$ から 0-1 の $n$ 次元バイナリからなるベクトル $e_1 \in \mathbb{R}^n$ を生成します。この $e_1$ を用いて Dropoutを適用した隠れ層 $h$ は

$$ \begin{align*} h = \sigma ( (x \bullet e_1) M_1 + b) \end{align*} $$

となります。 ここで $\bullet$ は要素積 $x \bullet y = \sum x_i y_i$ を表します。

同様に隠れ層 $h$ に対しても確率 $p_2$ で要素を0にします。すなわち 先と同様に $m$ 次元のベルヌーイ分布から $e_2 \in \mathbb{R}^m$ をサンプリングして一定の隠れ層ノードの値を0にします。すなわち

$$\hat{h} = h \bullet e_2$$

をノードの値であるとします。出力はこの値を用いて

$$\hat{y} = \hat{h} M_2$$

となります。
この $\hat{y}$ は $\bullet e = {\rm diag}(e)$ と変形できることを用いれば以下のように変形できます。

$$ \begin{align}\hat{y} & = (h \bullet e_2) M_2 \\&= ( \sigma \left( (x \bullet e_1) M_1 + b \right) \bullet e_2) M_2 \\&= ( \sigma (x ({\rm diag}(e_1) M_1) + b)({\rm diag}(e_2) M_2 ) \\&= \sigma \left( x\hat{W}_1 + b \right) \hat{W}_2 \end{align} $$

ここで  {\rm diag}(e_1) M_1 = \hat{W}_1, {\rm diag}(e_2) M_2 = \hat{W}_2 と定義しました。

以上を用いるとニューラルネットワークの出力は確率変数 $\hat{\omega} = \{ \hat{W_1}, \hat{W_2}, b\}$ を用いて

$$\hat{y} = f^{\hat{W_1}, \hat{W_2}, b}(x)$$

と記述できます. 


ながながとゴニョゴニョしましたが, 以上からdropout による mask のかけられた出力も重みの確率変数をもつネットワークの出力として表現できる, ということが確認できました。

Dropout の目的関数

これらの記号を用いてニューラルネットワークの目的関数を記述していきましょう。入力値と正解ラベルから誤差を計算するロス関数を $E$ とおき、ニューラルネットワークが最小化する真の関数を記述すると以下の用になります。

$$L_{dropout} = \frac{1}{M} \sum_{n \in N} E \left( f(x_n), y_n \right) + \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b ||$$

ここで $f(x_n)$ は $n$ 番目の入力データに対する出力値を表しています。また右辺第二項以降は重みに対する正則化 (weight decay)を表しています.

実際には $N$ 個すべてのデータを使うことは困難ですから、その中からある一定の大きさのサンプル $S$ を取得します。Dropout ではそれと同時にネットワークに対する mask をかけて出力を作ります。したがってネットワークの出力は確定値 $f(x)$ ではなく、確率変数 $\omega =\{ \hat{W_1}, \hat{W_2}, b\}$ によって決定する確率的な出力となります。よって Dropout を用いたニューラルネットワークのロス関数は

$$L_{dropout} = \frac{1}{M} \sum_{i \in S} E \left( f^{ \{ \hat{W_1^i}, \hat{W_2^i}, b \} }(x), y_i \right) + \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b ||$$

となります. ここで $\hat{W_1^i}, \hat{W_2^i}$ はサンプルされた $i$ のデータに対する dropout のマスクがかけられた重みを表しています。

回帰問題においては, 関数 $E$ は定数部分を除いて負の対数尤度関数で書き換えることができます。すなわち

$$E (f(x),y) = \frac{1}{2} \| y - f(x) \|^2 = - \frac{1}{\tau} \log p(y | f(x)) + {\rm const}$$

と変形できます。 ここで尤度関数は $p(y | f(x)) = N(y;f(x), \tau^{-1}I)$ のガウス分布で, $\tau$ は精度パラメータです.
今回は回帰問題を考えましたが、分類問題においても同様に負の対数尤度を用いて定式化することが可能です. (その場合 $\tau = 1$ となります.)

これよりロス関数は

$$L_{dropout} =  - \frac{1}{M} \sum_{i \in S} \frac{1}{\tau} \log p( y_i | f^{ \{ \hat{W_1^i}, \hat{W_2^i}, b \} }(x)) + \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b ||$$

つぎに確率変数 $\hat{\omega}$ について考えます。この集合はネットワークの重みという確定的な値とdropout による確率変数の部分をあわせたものでした。それを明示的に記述すると

$$\hat{\omega_i} = \{ \hat{W_1}, \hat{W_2}, b\} = \{{\rm diag}(e_1^i)M_1, {\rm diag}(e_2^i)M_2, b \} := g(\theta, \hat{e_i})$$

と書き換えることができます。ここで $\theta = \{M_1, M_2, b\}$ は確定的な値を要素に持つ集合と定義し、$\hat{e_i} = \{ e_1^i, e_2^i \}$ は $i$ 番目のミニバッチのサンプルによって作成される dropout のマスクを要素に持つ、確率変数の集合であると定義します。 また関数 $g$ は2つの集合 $\theta, \hat{e_i}$ から $\hat{\omega_i}$ をつくる射影であるとします。
また $1 \le i \le N$ にたいして $e_1^i \sim p(e_1)$, $e_2^i \sim p(e_2)$です. ここで, $p(e_j)\ (j=1,2)$ はそれぞれ確率 $p_j$ のベルヌーイ分布の積であるとします。
これらを用いると dropout のロス関数は確率変数 $\hat{e_i}$ を用いて

$$L_{dropout} = - \frac{1}{M\tau} \sum_{i \in S} \log p(y_i | f^{g(\theta, \hat{e_i})}(x_i)) + \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b ||$$

となります。この目的関数の重み $\theta := \{M_1, M_2, b\}$ に関する勾配は

$$\frac{\partial}{\partial \theta} L_{dropout} = - \frac{1}{M\tau} \sum_{i \in S} \frac{\partial}{\partial \theta} \log p(y_i | f^{g(\theta, \hat{e_i})}(x_i)) + \frac{\partial}{\partial \theta} \left( \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b || \right)$$

となり, この勾配を逆伝搬させてネットワークの重みを更新します.

ところでこれは先程の変分推論の式ととても良く似ています. 再掲すると

$$\hat{{\rm VI}}_{MC} = - \frac{N}{M} \sum_{i \in S} \log p(y_i | f^{g(\theta, \epsilon_i)}(x_i)) + {\rm KL} (q_{\theta}(w) | p(w))$$

であり, これは ${\rm KL} (q_{\theta}(w) | p(w)) = N\tau \left( \lambda_1 ||M_1|| + \lambda_2 || M_2 || + \lambda_3 || b || \right)$ とおけば

$$L_{dropout} = \frac{1}{M\tau} \hat{{\rm VI}}_{MC}$$

となりdropout と変分推論は厳密に一致します!!


厳密に一致するということは, dropout を用いて学習したネットワークは、変分推論による重みの事後分布となっているということです。

すなわち dropout を用いるのはベイズ学習をしていることにほかならない ということです。ですから学習によって得られている値は確定的な値ではなく、 dropout と組み合わせることで重みの事後分布を計算できます。

これによって。学習済みのネットワークにとあるdropoutをかけた出力は重みの事後分布から一つの重みをサンプルしているという風に解釈することが可能です。したがって、複数のdropoutの係数を用いて複数の出力を生成しそのばらつきを用いると, 事後分布の分散を推定することも可能です.

この枠組では, ニューラルネットワークの出力は確率過程からのサンプルにほかならないため、出力がどの程度信頼できるか(ばらつきをもっているか) を見積もることが可能になっています。すごい。

数値実験

簡単な人口データの回帰問題を解いてみます.

実験に用いたコードは以下においています。

github.com

条件

トレーニングデータは1次元の100個のデータで, ターゲット変数は $f(x) = x + \sin 5x$ に平均0, 分散 0.1 のガウスノイズを載せて作成します.

ネットワークは5層で各隠れ層が512次元としました。dropoutは確率 0.5 のベルヌーイ分布を用います。

relu

活性化関数を relu として 1000 epoch 計算したものの出力が下のグラフです。

f:id:dette:20170814172111p:plain

青い点がデータで, 橙色で事後分布からサンプルされたネットワークの出力を表しています。 これを見ると, データが存在する部分の分散は小さくなっていますがデータがない部分に行くに連れて分散が大きくなっていることがわかります. これは事後分布の分散が大きくなっているためです。このように従来のニューラルネットの予測値に加えて、その信頼度が出せていることが確認できます。

tanh

次に活性化関数を tanh に変えて実験したものが以下の図です。

f:id:dette:20170814174525p:plain

tanh に関してはreluのときのようにデータがない部分で極端に分散が大きくなりませんでした。またデータへの当てはまりも tanh の方は積極的に行っておらず、$x$ 軸性の領域ではほとんど 0 を予測しておりあたかも「データが足りないよ」と言っているようです。 

これは relu のほうが勾配消失に強く学習が素早く進むため、同じ epoch 数でもデータへ強くフィッティングしていくことが原因と考えられます。

また範囲外への予測も relu のようにそれまでの傾きを踏まえた予測ではなく、だんだんと 0 に向かうように予測しているように見受けられますが、これも活性化関数の形の影響 (reluは非ゼロのとき単調増加関数となる) を受けていると考えられます。

ちなみに tanh でも epoch さえ増やせばデータへしっかり当てはめてくれるようで、 epoch 数を更に増やして 4000 回としたものが以下になります。

f:id:dette:20170814175210p:plain

これを見るとデータがある部分に関しては平均値はほぼデータ通りで、分散もこちらが設定した値 0.1に近い値となっていることが伺えます。 1000 epoch では epoch が足りていなかったみたいですね。

感想

今まで確定値しか扱えなかったところに、Dropoutという既存の手法とベイズ推定の枠組みを結びつけて行くところは面白いなと思って読んでいました。CNNやRNNへ応用すれば「この画像はよくわかんないけど猫」とか「絶対犬です!」みたいな信頼度を同時にだすネットワークが自然に作れそうで今後発展すると面白そうだなと思います。

おまけ - 動画バージョン

tanh

 

f:id:dette:20170814180735g:plain

うにょうにょしながら学習していくの面白い。