①非線形分類問題（２クラス分類） - 社畜エンジニア発掘戦線

第一問

設問１．「原点を中心とした半径4の円」と「原点」にガウシアンノイズを付加したデータセット作成せよ
設問２．活性化関数にソフトマックス、損失関数にクロスエントロピーを用いて出力y1、y2と損失量Eを求めよ
設問３．勾配降下法でパラメータを600回更新し、その出力をデータセットと同じグラフ上にプロットせよ
まとめ

ニューラルネットワークの基本問題も終盤にかかってきました。これまでの分類問題は、基本的に直線で境界線を引いて分類できるデータセットの組み合わせでした。しかし、世の中には「紆余曲折」という言葉があるように曲がりくねったデータセットを取り扱う場合がほとんどです。この場合、どのように曲がった境界線を作るようなニューラルネットワークを構築するのか、設問を通じて取り組んでいこうと思います。

設問１．「原点を中心とした半径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(0~2π)を使って書いてみましょう。

$x=sin(s) \quad y=cos(s)$

この媒介変数にガウシアンノイズをのせます。

s = np.random.uniform(0,2*np.pi,num_of_sam)
x1 = np.sin(s)*radius
x2 = np.cos(s)*radius

しかし、このままではきれいな円上にまばらな間隔でデータが並ぶだけです、こんな感じ。

この円に幅を持たせる必要があるので、そのノイズも設定してあげなければいけません。そこで再度ガウシアンノイズを使って幅を広げるためのノイズ変数を定義します。

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

0.9~1.1はいい感じの幅になる値です、だから中のデータセットに干渉しない程度ならいくらでもいいんですけどね。最後にこのx1とx2を結合します。

X_circle = np.c_[x1,x2]

これで中心データセット「X_center」と環状データセット「X_circle」が完成しました。ニューラルネットワークに突っ込むためにはこれらをひとつのデータセットにまとめる必要があるので、最後にXとしてこの2つのデータセットを縦方向に統合します。

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))

これはワン・ホットで線形分類問題のときに散々やりましたね。

ではグラフ化してみます。

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:20190304075453p:plain
おお〜、これこれ、いい感じです。長くなりましたが設問1は終了です。

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

さて、今回は「非線形」問題です。非線形と名の付くものに隠れ層あり、そのことを念頭に置いてニューラルネットワークを設計します。

隠れ層Hを1層追加しました。そのニューロンは今回は3つにしましょう。隠れ層の活性化関数にはシグモイド、分類問題なのでこれまで通り出力層にはソフトマックス関数を使います。損失関数もこれまで通りクロスエントロピーを用います。いつかの問題で書いたように、ワン・ホットの分類問題で出力層に使う活性化関数はソフトマックスが良いということでしたが、隠れ層の活性化関数はシグモイドじゃなくても良さげです（最近はReLU関数とか流行っているらしい）。ずっとシグモイド関数しか使ってなかったので、ここではシグモイドを使いますが。

#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)))

次にパラメータWとBの設定です。各層でのパラメータの型を確認してみます。
確認のために数式で書き下してみるとこんな感じ、ただしここでは行列を書きやすいように転置してるのでご注意。

ではそれぞれのパラメータをコードに落とします。

#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)

ここは記述そのままですね。ニューラルネットワークがやや複雑なのでXからそれぞれの型を確認しておきます。

print(X.shape) #(80, 2)
print(H.shape) #(80, 3)
print(Y.shape) #(80, 2)
print(E.shape) #()

サンプル数は40で入力Xはx1とx2の2つなので#(80, 2)、隠れ層Hはニューロンがひとつ増えて#(80, 3)、出力Yは2つなので#(80, 2)、損失関数はスカラーなので型なし、オッケーです。

設問2は終了です。

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

まずはパラメータの更新式を求めるわけですが、今回の問題には隠れ層があるのでもう一度微分式から確認してみます。そう言えば隠れ層は非線形回帰問題の1回しか扱っていないから全然慣れてない。ですが、流れはこれまでの計算と同じです。適時こちらの計算を利用してるのでご参照。
活性化関数とクロスエントロピーの微分 - 社畜エンジニア発掘戦線

まずは出力層のパラメータ微分から、WとBを一気にやってしまいます。

${\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) }$

ここまでは前回までよくやっていたのでいいんですが、次の隠れ層の逆伝搬は慎重にいきます。まずはW1だけで見てみます。

${\displaystyle \frac{\partial E}{\partial W_1} = \frac{\partial H}{\partial W_1} \frac{\partial Y}{\partial H} \frac{\partial E}{\partial Y} }$

2層戻るので連鎖率は3段階、次に後ろの2つだけを計算します。ここもソフトマックスとクロスエントロピーの計算なのでΔが使えます。

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

次に隠れ層の微分計算です。ここにはシグモイド関数がかかっているので注意です。よくこの計算を忘れてプログラムの結果が合わずにめっちゃ悩むことあるんですよね…。

${\displaystyle \frac{\partial H}{\partial W_1} = \frac{\partial }{\partial W_1} sigmoid(XW_1+B_1) = X \cdot H(1-H) }$

最後に全部を統合して微分式を求めます。

${\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) }$

なんかめっちゃ計算しやすくなった気がする、やっぱΔすげぇよ。合わせてB1についても計算します。

${\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) }$

おお、できた。最後にこの微分式を使って更新式も求めておきます。と言っても勾配降下法なので今までと同じですが。ただし、出力層と隠れ層2つあるので合計4つのパラメータ更新が必要です（添字が見にくいのはご勘弁）。

${\displaystyle W_{1-new} = W_{1-old} - \frac{\partial E}{\partial W_1} \qquad B_{1-new} = B_{1-old} - \frac{\partial E}{\partial B_1} }$

${\displaystyle W_{2-new} = W_{2-old} - \frac{\partial E}{\partial W_2} \qquad B_{2-new} = B_{1-old} - \frac{\partial E}{\partial B_2} }$

ではコードに落とします。イテレーション回数は600回、学習率は0.008とします（毎度、何回か動かしてみてしっくりきそうな値をチョイスしています）。

#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)

learning_rate = 0.008
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
    W1 = W1 - learning_rate*dW1
    B1 = B1 - learning_rate*dB1
    W2 = W2 - learning_rate*dW2
    B2 = B2 - learning_rate*dB2

パラメータの更新が4つになったので長くなっていますが、構成は今までと同じです。

ではグラフ化してみます。
まずはグリッドに切るところから。ここも今までと同じです。

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]

グラフを描画します。

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:20190304080020p:plain
おお、うまくいきました。

めでたく設問3が終了です。

まとめ

いつもはここで出力だけをGIFアニメ化するんですが、確認したいこともあるので損失量のグラフと一緒にアニメ化してみます。
f:id:sutokun:20190304080151g:plain f:id:sutokun:20190304080214g:plain
何度かトライしてみると、「うまくいくとき」と「いかないとき」があります。体感では30％ぐらい失敗かな。

うまくいかなかったとき。
f:id:sutokun:20190304080918g:plain f:id:sutokun:20190304080943g:plain
上のグラフほど損失量が収束していません、最後の方にちょっと良くなり始めてる感じはありますが。

同じコードのプログラムですが損失量の収束に差があります。これは損失関数上の「局所解」にハマっていると思われます。（ちなみに学習の出発点はWをランダムな行列で定義しているので毎回変わります）以前に線形回帰問題で損失量の可視化にトライしました。その問題では、損失量が二次関数で表現された損失関数上を転がって最小値へ向かう様子が確認できました。
→　③線形回帰問題（損失関数の可視化） - 社畜エンジニア発掘戦線

線形回帰問題では設定が簡単だったので「極値＝最小値」となっていました。が、今回のような問題では損失関数がどのような形状をしているのかは分かりません。勾配降下法において損失量はイテレーションごとに、その位置に対して、その位置よりも低い方向へ向かいます。今回の損失関数がいくつかの局所解を持っている場合、その極値にハマってしまっている可能性が高いです。