ニューラルネットワークを実装する [Part 4 バックプロパゲーションのベクトル化]

はじめに

前回に引き続き、Peters' NoteのPart 4 Vectorizationを見ていきます。

今回の記事では2クラス分類問題に対して、

バックプロパゲーションをベクトルで扱う方法
Gradient Checking
モーメンタム法

を順に実装していきます。

はじめに
目次
今回扱うニューラルネットワーク
使用するデータセットの用意
ベクトル化
- 順伝播
  - 隠れ層の出力
  - 出力層の出力
- 逆伝播
  - 出力層での逆伝播
  - 隠れ層での逆伝播
Gradient Checking
モーメンタム法を使った最適化
最急降下法を使ったらどうなるの？
GitHub
関連記事

今回扱うニューラルネットワーク

今回は、下図のような入力層、1層の隠れ層、出力層がつながったモデルを扱います。2次元の入力を受け取り、3次元の隠れ層、2次元の出力層へとつながっていきます。
前回までと異なり、今回はバイアスも定義されています。

f:id:kanohk:20171104154144p:plain — 出典: https://peterroelants.github.io/posts/neural_network_implementation_part04/#Vectorization

出力層の活性化関数にsoftmax関数を使うと、このモデルは2クラス分類問題に使用することができます。

使用するデータセットの用意

いつも通り使用するデータセットを生成します。今回は、下図のような、青クラスの外側を赤クラスが取り囲んだようなデータセットを使用します。
f:id:kanohk:20171104154523p:plain

# データセットを用意
X, t = sklearn.datasets.make_circles(n_samples=100, shuffle=False, factor=0.3, noise=0.1)
T = np.zeros((100, 2))
T[t == 1, 1] = 1
T[t == 0, 0] = 1

# 可視化のために2クラスを別の変数に分ける
x_red = X[t == 0]
x_blue = X[t == 1]

print('shape of X: {}'.format(X.shape)) # shape of X: (100, 2)
print('shape of T: {}'.format(T.shape)) # shape of T: (100, 2)

scikit-learnのmake_circlesという関数を使うと、今回使いたいデータを簡単に用意することができます。
x_redとx_blueをプロットすると先ほど示した図が得られます。

ベクトル化

順伝播

隠れ層の出力

モデルへの入力は2次元のベクトル $\mathbf{X}$ で与えられます。隠れ層の値 $\mathbf{H}$ は3次元のベクトルで、重み $\mathbf{W}_h$ とバイアス $\mathbf{b}_h$ を使うと次のように表すことができます。
$\displaystyle \mathbf{H} = \sigma (\mathbf{X} \cdot \mathbf{W}_h + \mathbf{b}_h)$

ただし、 $\sigma$ はロジスティック関数、 $\mathbf{W}_h$ と $\mathbf{b}_h$ は
$\displaystyle \begin{align} \mathbf{W}_h = \begin{bmatrix} w_{h11} & w_{h12} & w_{h13} \\ w_{h21} & w_{h22} & w_{h23} \end{bmatrix} && \mathbf{b}_h = \begin{bmatrix} b_{h1} & b_{h2} & b_{h3} \end{bmatrix} \end{align}$
で定義されます。

この数式をコードに落とすとこのようになります。

# ロジスティック関数
def logistic(z):
    return 1 / (1 + np.exp(-z))

# 隠れ層の出力
def hidden_activations(X, Wh, bh):
    return logistic(X.dot(Wh) + bh)

出力層の出力

上記と同様に、出力層の値 $\mathbf{Y}$ も、重み $\mathbf{W}_o$ とバイアス $\mathbf{b}_o$ 、Softmax関数 $\varsigma$ を使用して次のように記述できます。
$\displaystyle \mathbf{Y} = \varsigma (\mathbf{H} \cdot \mathbf{W}_o + \mathbf{b}_o)$

# softmax
def softmax(z):
    return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)

# 出力層の出力
def output_activations(H, Wo, bo):
    return softmax(H.dot(Wo) + bo)

以上をもとに、ニューラルネットワーク全体の関数と、入力からクラスを予測する関数を定義することができます。

# ニューラルネットを定義
def nn(X, Wh, bh, Wo, bo):
    return output_activations(hidden_activations(X, Wh, bh), Wo, bo)

# 入力からクラスを予測する関数
def nn_predict(X, Wh, bh, Wo, bo):
    return np.around(nn(X, Wh, bh, Wo, bo))

逆伝播

出力層での逆伝播

式が長くなるので、細かい計算は元のノートを見ていただきたいのですが、コスト関数にクロスエントロピー $\xi (\mathbf{T}, \mathbf{Y})$ を使用した時に、バックプロパゲーションで使用するコスト関数の微分の値は次のようになります。
$\displaystyle \mathbf{\delta}_{o} = \frac{\partial \mathbf{\xi}}{\partial \mathbf{Z}_o} = \mathbf{Y} - \mathbf{T}$
ただし、 $\mathbf{Z}_o$ は出力層(ソフトマックス層)への入力です。

コスト関数と $\mathbf{\delta}_o$ (Python内ではerror_output関数)は次のように定義できます。

# コスト関数
def cost(Y, T):
    return -np.multiply(T, np.log(Y)).sum()

# 出力層の誤差
def error_output(Y, T):
    return Y - T

さらに、コスト関数の重みでの微分は、(こちらも計算過程は一切省略してしまいますが。。）次のように書けます。
$\displaystyle \frac{\partial \xi}{\partial W_o} = H^T \cdot (Y-T) = H^T \cdot \delta_{o}$

この値はヤコビアン行列なので、 $J_{W_o}$ と簡単に記述します。

同様にバイアスに関してもヤコビアンが次のように記述できます。
$\displaystyle J_{b_o} = \frac{\partial \xi}{\partial \mathbf{b}_o} = \sum_{i=1}^N \delta_{oi}$

この二つのヤコビアンを計算する関数は、次のように実装できます。

# コスト関数のWoでの微分
def gradient_weight_out(H, Eo):
    return H.T.dot(Eo)

# コスト関数のboでの微分
def gradient_bias_out(Eo):
    return np.sum(Eo, axis=0, keepdims=True)

隠れ層での逆伝播

再び計算式はノートに譲るとして、計算結果だけを紹介すると、ヤコビアンは次のように記述できます。
$\displaystyle \delta_{h} = \frac{\partial \xi}{\partial Z_h} = H \circ (1 - H) \circ [\delta_{o} \cdot W_o^T$ ]
$\displaystyle J_{W_h} \frac{\partial \xi}{\partial W_h} = X^T \cdot \delta_{h}$
$\displaystyle J_{b_h} = \frac{\partial \xi}{\partial \mathbf{b}_{h}} = \sum_{j=1}^N \delta_{hj}$
ただし、 $\circ$ は要素ごとの積を表しています。

途中の計算こそややこしいのですが、計算結果はシンプルで、簡単なコードに落とすことができます。

# 隠れ層の誤差
def error_hidden(H, Wo, Eo):
    # H * (1 - H) * (E, Wo^T)
    return np.multiply(np.multiply(H, (1 - H)), Eo.dot(Wo.T))

# コスト関数のWhでの微分
def gradient_weight_hidden(X, Eh):
    return X.T.dot(Eh)

# コスト関数のbhでの微分
def gradient_bias_hidden(Eh):
    return np.sum(Eh, axis=0, keepdims=True)

Gradient Checking

バックプロパゲーションの実装はバグを埋め込みやすく、実装が正しいかをチェックすることは非常に大事だそうです。
バックプロパゲーションで計算した微分の値と、数値微分の値が大きく異なっていないかを調べることで、実装の正しさをチェックできます。

数値微分とは、次の式で微分の値を近似することをいいます。
$\displaystyle \frac{\partial \xi}{\partial \theta_i} = \frac{f(X, \theta_1, \cdots, \theta_i+\epsilon, \cdots, \theta_m) - f(X, \theta_1, \cdots, \theta_i-\epsilon, \cdots, \theta_m)}{2\epsilon}$
ただし、 $\epsilon$ は十分に小さな値でなければいけません。

では実際にバックプロパゲーションで計算した微分と、数値微分の値が近いものかを調べてみます。
まずはバックプロパゲーションで微分を計算します。順方向に出力を計算した後、微分を値を逆方向に計算しています。

# 重みとバイアスの初期化
init_var = 1
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var

# バックプロパゲーションで微分を計算
# (1) 順方向に出力を計算
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)

# (2) 逆方向に微分を計算
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)

Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

次に各パラメータに対して数値微分を計算して、バックプロパゲーションでの微分の値と近い値かをnumpyのisclose関数を使用して調べます。

# パラメータを1つのリストにまとめておく
params = [Wh, bh, Wo, bo]

# バックプロパゲーションで計算した微分
grad_params = [JWh, Jbh, JWo, Jbo]

# 数値微分に使うイプシロン
eps = 0.0001

# Gradient checking
for p_idx in range(len(params)):
    for row in range(params[p_idx].shape[0]):
        for col in range(params[p_idx].shape[1]):
            # 対象のパラメータ行列をコピー
            p_matrix_min = params[p_idx].copy()
            p_matrix_min[row, col] -= eps
            p_matrix_plus = params[p_idx].copy()
            p_matrix_plus[row, col] += eps
            
            # paramsをコピーし、更新した行列で置き換える
            params_min = params[:]
            params_min[p_idx] = p_matrix_min
            params_plus = params[:]
            params_plus[p_idx] = p_matrix_plus
            
            # 数値微分を計算
            grad_num = (cost(nn(X, *params_plus), T) - cost(nn(X, *params_min), T)) / (2 * eps)
            
            # np.iscloseでバックプロパゲーションでの微分と数値微分の値が近いか確認
            if not np.isclose(grad_num, grad_params[p_idx][row, col]):
                raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(
                    float(grad_num), float(grad_params[p_idx][row, col])))

print('No gradient errors found')

このコードを実行すると、No gradient errors foundと表示され、全てのパラメータについて微分の値が近いことが確認できました。

モーメンタム法を使った最適化

以前は最急降下法を使用してコスト関数の最適化を行っていましたが、今回はモーメンタム法を使用します。後に見るように、最急降下法だと局所最適解にたどり着いてしまうことが多く、良い結果が得られないからです。

モーメンタム法では、次の式でパラメータを更新します。
$\displaystyle \begin{split} V(i+1) & = \lambda V(i) - \mu \frac{\partial \xi}{\partial \theta(i)} \\ \theta(i+1) & = \theta(i) + V(i+1) \end{split}$
$V(i)$ はイテレーション $i$ での速度で、その速度を使用してパラメータ $\theta$ を更新していきます。

# バックプロパゲーションで微分を計算する関数
# やっていることは先ほど示したコードと同じです
def backprop_gradients(X, T, Wh, bh, Wo, bo):
    H = hidden_activations(X, Wh, bh)
    Y = output_activations(H, Wo, bo)
    
    Eo = error_output(Y, T)
    JWo = gradient_weight_out(H, Eo)
    Jbo = gradient_bias_out(Eo)
    
    Eh = error_hidden(H, Wo, Eo)
    JWh = gradient_weight_hidden(X, Eh)
    Jbh = gradient_bias_hidden(Eh)
    
    return [JWh, Jbh, JWo, Jbo]

# 速度を更新する関数
def update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate):
    # ls_of_params = [Wh, bh, Wo, bo]
    # Js = [JWh, Jbh, JWo, Jbo]
    Js = backprop_gradients(X, T, *ls_of_params)
    return [momentum_term * V - learning_rate * J for V, J in zip(Vs, Js)]

# パラメータを更新する関数
def update_params(ls_of_params, Vs):
    # ls_of_params = [Wh, bh, Wo, bo]
    # Vs = [VWh, Vbh, VWo, Vbo]
    return [P + V for P, V in zip(ls_of_params, Vs)]

いよいよバックプロパゲーションを使ってネットワークの最適化を行っていきます。
update_params関数を使ったパラメータの更新を、300イテレーション回します。

# バックプロパゲーションの実行

# 重みとバイアスの初期化
init_var = 0.1
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var

learning_rate = 0.02
momentum_term = 0.9

Vs = [np.zeros_like(M) for M in [Wh, bh, Wo, bo]]

# イテレーション部分
num_of_iterations = 300
lr_update = learning_rate / num_of_iterations
ls_costs = [cost(nn(X, Wh, bh, Wo, bo), T)] # イテレーションごとのコストを保存しておく

for i in range(num_of_iterations):
    Vs = update_velocity(X, T, [Wh, bh, Wo, bo], Vs, momentum_term, learning_rate)
    Wh, bh, Wo, bo = update_params([Wh, bh, Wo, bo], Vs)
    ls_costs.append(cost(nn(X, Wh, bh, Wo, bo), T))

ls_costsに保存したイテレーションごとのコスト関数の値を可視化すると次のようになります。最初は70程度だった値が、200回を超えたあたりから0に近い小さな値に収束していく様子が確認できます。
f:id:kanohk:20171104165858p:plain

学習したパラメータを使って、決定境界を可視化すると次のようになります。うまく入力を分類できている様子が確認できます。
f:id:kanohk:20171104170007p:plain

隠れ層は3次元でしたので、隠れ層の値も3次元空間上にプロットすることができます。それがこちらです。隠れ層の時点で3次元空間上では青クラスと赤クラスがうまく分類できているようです。
f:id:kanohk:20171104170111p:plain

最急降下法を使ったらどうなるの？

今回はモーメンタム法を使って最適化を行いましたが、前回までと同様の最急降下法を使うとどのような結果になるのでしょうか？
最急降下法では、 $W(i + 1) = W(i) - \Delta W(i)$ でパラメータの更新をするので、コードは次のようになるかと思います。

# モーメンタムでなく最急降下法を使用して最適化してみる
def gradient_decent_params_update(X, T, Wh, bh, Wo, bo, learning_rate):
    JWh, Jbh, JWo, Jbo = backprop_gradients(X, T, Wh, bh, Wo, bo)
    return [W - learning_rate * dW for W, dW in zip([Wh, bh, Wo, bo], [JWh, Jbh, JWo, Jbo])]

# バックプロパゲーションの実行
# 重みとバイアスの初期化
init_var = 0.1
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var

learning_rate = 0.2

# イテレーション部分
num_of_iterations = 300
lr_update = learning_rate / num_of_iterations
ls_costs = [cost(nn(X, Wh, bh, Wo, bo), T)] # イテレーションごとのコストを保存しておく

for i in range(num_of_iterations):
    learning_rate -= lr_update
    Wh, bh, Wo, bo = gradient_decent_params_update(X, T, Wh, bh, Wo, bo, learning_rate)
    ls_costs.append(cost(nn(X, Wh, bh, Wo, bo), T))

再び最適化の過程を可視化すると次のようになります。
f:id:kanohk:20171104170501p:plain

振動しながらも徐々にコスト関数は小さくなっていきますが、最終的に70-80程度までしかコストを小さくできていません。
Jupyter Notebookで何回か上記コードを実行して試してみましたが、何度やっても同様にある程度の値までコストが小さくなったらそれ以上は小さくなりませんでした。おそらく、局所最適解につかまってしまっているのだと思います。

念のため、決定境界も可視化してみましたが、やはり全然うまく分類できていない様子が確認できます。
f:id:kanohk:20171104170721p:plain

GitHub

可視化部分のコード等も含めたJupyter Notebookは下記においてあります。
notebooks/Part4.ipynb at master · kanoh-k/notebooks · GitHub

これまでの記事

今日も窓辺でプログラム

外資系企業勤めのエンジニアが勉強した内容をまとめておくブログ