社畜エンジニア発掘戦線

駆けだしAIエンジニア

②非線形分類問題（パラメータの更新）

第二問

前回は非線形分類問題ということで、ドーナツ型のデータセットを分類しました。しかし、勾配降下法を用いると損失関数の局所解にハマり学習が進まなくなる問題が発生しました。今回は勾配降下法以外のパラメータ更新方法でこの局所解にハマっちゃう問題がどう改善されるのかトライしていきます。

設問１．「原点を中心とした半径4の円」と「原点」にガウシアンノイズを付加したデータセット作成せよ

ここは前回と同じです。そのままコピペしてしまいましょう。まずはライブラリのインポートから。

import numpy as np
import matplotlib.pyplot as plt

#dataset
num_of_sam = 40
std_dv = 0.6
radius = 4

X_center = np.random.randn(num_of_sam,2)*std_dv

s = np.random.uniform(0,2*np.pi,num_of_sam)
noise = np.random.uniform(0.9, 1.1, num_of_sam)
x1 = np.sin(s)*radius*noise
x2 = np.cos(s)*radius*noise
X_circle = np.c_[x1,x2]

X = np.vstack((X_center,X_circle))

t_group1 = np.tile([0,1],(num_of_sam,1))
t_group2 = np.tile([1,0],(num_of_sam,1))
T = np.vstack((t_group1, t_group2))

そしてグラフ化、
f:id:sutokun:20190317011751p:plain
そらコピペなんでいい感じです、次へ進みましょう。

設問２．活性化関数にソフトマックス、損失関数にクロスエントロピーを用いて出力y1、y2と損失量Eを求めよ

ここも前回と同じです。今回はパラメータの更新手法にフォーカスしているのでフォワードプロップ（順伝搬計算）はすべて前回と同じです。

まずは関数の書き下し。

#function
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

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

def loss(y, t):
    return -np.sum(np.multiply(t, np.log(y)) + np.multiply((1 - t), np.log(1 - y)))

次に初期値を設定して、順計算を進めます。

#initial setting
W1 = np.random.randn(2,3)
B1 = np.random.randn(1,3)
W2 = np.random.randn(3,2)
B2 = np.random.randn(1,2)

H = sigmoid(np.dot(X,W1)+B1)
Y = softmax(np.dot(H,W2)+B2)
E = loss(Y, T)

こんな感じでこの設問は終了です。

設問３．損失関数を各層のパラメータで微分した微分式を求めよ

前の問題では損失関数まで出せたので、これを各パラメータで微分します。これも前回やったので結果だけ。

まずは数式から。

{\displaystyle \frac{\partial E}{\partial W_2} = \frac{\partial Y}{\partial W_2} \frac{\partial E}{\partial Y} = \frac{\partial }{\partial W_2} (H W_2 + B_2) \cdot \Delta = H \cdot (Y-T) }

{\displaystyle \frac{\partial E}{\partial B_2} = \frac{\partial Y}{\partial B_2} \frac{\partial E}{\partial Y} = \frac{\partial }{\partial B_2} (H W_2 + B_2) \cdot \Delta = 1 \cdot (Y-T) }

{\displaystyle \frac{\partial E}{\partial W_1} = \frac{\partial H}{\partial W_1} \frac{\partial Y}{\partial H} \frac{\partial E}{\partial Y} = X \cdot H(1-H) \cdot W_2 \cdot (Y-T) }

{\displaystyle \frac{\partial E}{\partial B_1} = \frac{\partial H}{\partial B_1} \frac{\partial Y}{\partial H} \frac{\partial E}{\partial Y} = 1 \cdot H(1-H) \cdot W_2 \cdot (Y-T) }

コードに落とします。

#back propagation
dW2 = np.dot(H.T,Y-T)
dB2 = np.sum(Y-T, axis=0, keepdims=True)
dW1 = np.dot(X.T,H*(1-H)*np.dot(Y-T,W2.T))
dB1 = np.sum(H*(1-H)*np.dot(Y-T,W2.T), axis=0, keepdims=True)

以上で設問は終了です。

設問４．モーメンタム法でパラメータを600回更新し、その出力をデータセットと同じグラフ上にプロットせよ

新キャラが出てきました。モーメンタム法とは勾配降下法とは異なるパラメータ更新手法です。考え方としては、勾配降下法が「今の位置に対して低い方へ向かう」だったのに対して、モーメンタム法は「重力？をもって低い方へ落ちる」です。例えば、小さな谷を持つ損失関数上を損失量が動くとき、勾配降下法だとその次の小さな山を超える前に更新が止まってしまうと、その位置よりも低いのは「小さな谷」側なので、その結果として損失量は「小さな谷」にハマってしまいます。これがモーメンタム法だと損失量の動きに重力のような力が働くので、仮に小さな山を超える前に更新が止まってしまっても、その慣性で山を超えることができます。
f:id:sutokun:20190309125345j:plain
※勾配降下法はstochastic gradient descentらしいので略してSGDと書くことが多いです

ではどのように重力を付加させるのか数式を確認してみます。モーメンタム法では新たに「速度（Verocity）」というパラメータを追加します。速度を更新して、パラメータを更新する、という2段階アップデートです。この速度には学習率μと合わせて摩擦に相当する減衰率γが必要になります。

まずは速度の更新式、Pはパラメータの略でWとかBが入ります。微分式はこれまでと同じ損失関数を各層のパラメータで微分した式になります。

{\displaystyle V_{new} = \gamma V_{old}-\mu \frac{\partial E}{\partial P} }

続いてパラメータの更新式です。

{\displaystyle P_{new} = P_{old}+V_{new} }

では、前の設問で出した微分式とこの更新式を合わせてコードに書いていきます。まずはモーメンタムの速度を定義します。各パラメータの行列の型と同じにして、最初はすべてゼロでいきましょう。

velocity_W1 = np.zeros_like(W1)
velocity_B1 = np.zeros_like(B1)
velocity_W2 = np.zeros_like(W2)
velocity_B2 = np.zeros_like(B2)

この速度パラメータに加えて減衰係数γも必要です。これらをすべて含めて初期値を下記のように書きます。なお、減衰率は一般的に約0.9ぐらいの値を設定するそうです。理由はよく分からんですが、そうみたいです。

#initial parameter setting
W1 = np.random.randn(2,3) #重みとバイアス
B1 = np.random.randn(1,3)
W2 = np.random.randn(3,2)
B2 = np.random.randn(1,2)

velocity_W1 = np.zeros_like(W1) #速度パラメータ、全部ゼロ
velocity_B1 = np.zeros_like(B1)
velocity_W2 = np.zeros_like(W2)
velocity_B2 = np.zeros_like(B2)

learning_rate = 0.008 #学習率
momentum_term = 0.9 #減衰率
E_save = [] #損失量のセーブ

では実際にイテレーションを回すコードを書いていきます。上で計算した更新式に沿ってコードに落としていきます。大まかな流れはこれまでと同じです。

#iteration
num_of_itr = 600
for i in range(num_of_itr):
    #forward propagation
    H = sigmoid(np.dot(X,W1)+B1)
    Y = softmax(np.dot(H,W2)+B2)
    E = loss(Y, T)
    E_save = np.append(E_save, E)

    #back propagation
    dW2 = np.dot(H.T,Y-T)
    dB2 = np.sum(Y-T, axis=0, keepdims=True)
    dW1 = np.dot(X.T,H*(1-H)*np.dot(Y-T,W2.T))
    dB1 = np.sum(H*(1-H)*np.dot(Y-T,W2.T), axis=0, keepdims=True)

    #update parameter
    velocity_W1 = momentum_term*velocity_W1-learning_rate*dW1
    W1 = W1+velocity_W1

    velocity_B1 = momentum_term*velocity_B1-learning_rate*dB1
    B1 = B1+velocity_B1

    velocity_W2 = momentum_term*velocity_W2-learning_rate*dW2
    W2 = W2+velocity_W2

    velocity_B2 = momentum_term*velocity_B2-learning_rate*dB2
    B2 = B2+velocity_B2

ではグラフ化していきます。まずはグリッドに切って、プロット、一気にいきましょう。

#grid setting
grid_range = 10
resolution = 50
x1_grid = x2_grid = np.linspace(-grid_range, grid_range, resolution)

xx, yy = np.meshgrid(x1_grid, x2_grid)
X_grid = np.c_[xx.ravel(), yy.ravel()]

H_grid = sigmoid(np.dot(X_grid, W1)+B1)
Y_grid = softmax(np.dot(H_grid, W2)+B2)
Y_predict = np.around(Y_grid)

out_connect = np.hstack((X_grid,Y_predict))
blue_group = out_connect[out_connect[:,2]==1]
red_group = out_connect[out_connect[:,3]==1]

#show graph
plt.figure()
plt.plot(blue_group[:,0],blue_group[:,1],'o',alpha=0.3,color='blue')
plt.plot(red_group[:,0],red_group[:,1],'o',alpha=0.3,color='red')

plt.plot(X_center[:,0],X_center[:,1], 'o',color='red')
plt.plot(X_circle[:,0],X_circle[:,1], 'o',color='blue')

plt.show()

f:id:sutokun:20190317011904p:plain

おお、すばらしい、分類はできてそうです。
無事に設問4を終了できました。

まとめ

今回は勾配降下法以外のパラメータ更新方法でどのような違いが見られるのかを確認しました。設問ではモーメンタム法の収束の様子しかなかったので勾配降下法と比べてどれぐらい早いのかよく分かりません。下に同じデータセットで降下法とモーメンタム法の違いを実際の分類の様子と損失関数のグラフでアニメ化しています。左が降下法で右がモーメンタム法です。損失量のグラフは青と緑がそれぞれ降下法とモーメンタム法です。


f:id:sutokun:20190317012812g:plainf:id:sutokun:20190317012832g:plain
f:id:sutokun:20190317012859g:plain

あきらかにモーメンタム法の方が降下法よりも早く収束していることが分かります。おそらく、モーメンタム法はうまく局所解を抜け出して最小値へたどり着いているんだと思います。てかこんなに違いが見えるんですね、へーすごい。


sin関数をフィットさせた非線形回帰問題をトライしたとき、隠れ層では多数のシグモイド関数が足し合わされる様子を確認しました（あー、隠れ層の可視化までしてなかった、忘れてた、まあいいか、またやります）。
→①非線形回帰問題（データフィッティング） - 社畜エンジニア発掘戦線

分類問題では隠れ層でどのような動きが繰り広げられているのか気になるところです。その辺を下記にまとめているのでまた参照ください。あとコーディングの改良もしています。
→隠れ層の可視化 - 社畜エンジニア発掘戦線



これでディープラーニングの基礎問題は終了です。次回以降は少し内容を発展させて画像認識（MNISTの手書き数字）問題にトライしたいと思います。ただ、MNISTのデータセットを読み込むのにどうにもPYTHONISTAだとうまくいかず、次回はPC(OSはMac)でのトライになります。PYTHONISTAはお役御免や、すまんな。

なんやかんやここまでPYTHONISTAを使って進めることができました。処理落ちするとか文句を言いつつもムチを打てばアヒぃと言いながら頑張ってくれたPYTHONISTAには感謝です。

次の問題からはいったん仕切り直してこちらのページにまとめていきたいと思います。
週末のディープラーニング - 社畜エンジニア発掘戦線


全体コード
github.com

元の記事
PYTHONISTA3を使って機械学習（ディープラーニング） - 社畜エンジニア発掘戦線

Twitter
世界の社畜 (@sekai_syachiku) | Twitter