社畜エンジニア発掘戦線

駆けだしAIエンジニア

週末のディープラーニング

プログラミング

おつかれさまです。

これまでPYTHONISTAを使ってディープラーニングの基礎問題をトライしてきました。これらの問題を通してディープラーニングの「デ」ぐらいは理解できたんじゃないかな〜と自負している今日このごろです（すいません、まだまだ初心者です）。
shachicode.hatenablog.com

もっと記事にしたい小ネタはあるんですが、早く発展的な問題もトライしていきたいってことで、いったん仕切り直しました。

CONTENTS

開発環境

ここからはPYTHONISTAじゃインストールできないライブラリや、iPadの計算能力の問題もあり、PC（OSはMac）を使って発展問題にトライしていきたいと思います。

使っているMacのスペックはこんな感じです、「このマックについて」ってところをポチー。

プロセッサ：2.4GHz Inter Core i5
メモリ：4GB 1600MHz DDR3
グラフィックス：Inter Iris 1536MB

なにせディープラーニングは初学者なため、このスペックが良いのか悪いのか分かりません。ただ、Webの記事とかを読んでると明らかに悪そう。「グラフィックボードはウンGBのやつを積んでだなぁ！！」っていっぱい書いてて「あぁ〜（哀愁」となっています。

とにかく、できるところまでこのPCを使ってみて、時間がかかってしょうがない！ってなったらその時に考えることにします。まぁ昼は仕事に行ってるし、毎日夜は寝てるし、その間に頑張ってくれればなんとかなるんちゃう？（楽観

ただ、今アメリカに住んでるのでグラボの調達とかできるんかな…？しかもMacのラップトップだからグラボ増設とかそもそも無理じゃね？いやいや、その時に考えよう。


まだまだディープラーニングのアルゴリズムにはnumpyだけでチャレンジしようと思っているんですが、いずれどこかのタイミングでフレームワークを導入することになるかと思っています。その日に備えてGoogle先生の「tensorflow」をインストールしました。しかし、tensorflowを使うにはそれなりのPCスペックを要求されるらしく、たぶん使いこなすのは無理かな〜とも身構えています。

その他、エディターとかも更新しました。

↓こちらにまとめています。
開発環境のアップデート - 社畜エンジニア発掘戦線


（追記．6/4）
ニューラルネットワークの構築を進めていると、マジで計算時間がかかってきました。PCのスペックをどうするか検討中です。。。

練習問題
ディープラーニングのコツなど
Pythonのコードについてメモ


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

隠れ層の可視化

おつかれさまです。

非線形分類問題では隠れ層を用いて直線では分類できないような2クラスの分類をトライしました。
→②非線形分類問題（パラメータの更新） - 社畜エンジニア発掘戦線

出力層で確認できる結果としては、確かに2クラスの分類ができていましたが、「じゃあ隠れ層ではどうなっとんねん！」っていうのが素朴な疑問だと思います。

この問題では隠れ層（H）を3つのニューロンで構成しました、つまりH1、H2、H3の3次元グラフで可視化することができます。今回はこの隠れ層でどんなふうに変化が起きているのか確認してみたいと思います。使うコードは上記のリンクの問題の完成コードから引用。ちなみに3Dグラフを使うのでPYTHONISTAでは確認できません、PC（Mac）を使ってトライしてみます。

f:id:sutokun:20190318124520j:plain:w500

合わせて、ちょろっとコードの改善もトライします。基本的に設問の中ではできるだけ数式に忠実にコード化してたんですが、繰り返しもいちいちコード化してブサイクなところも多々あり（まぁ分かりやすいんだけど）、そこらへんも一度まとめてきれいにしたいと思います。

ライブラリのインポートからデータセットまで

まずはライブラリのインポート、今回は3Dグラフを使うのでそれ用の呪文を唱えます。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D #3D呪文

次にデータセット、ここはそのまま使います。まとめれなくもないんですが可読性が下がりすぎるのでやめました。

#data set
num_of_sam = 40
radius = 4
std_dv = 0.6

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))
関数の定義1

それから関数の定義です、シグモイド、ソフトマックス、クロスエントロピーをそれぞれコードに落としておきます。

#function1
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)))
関数の定義2

いつもはここから初期値を設定して、ニューラルネットワークの計算を進めていたんですが、このニューラルネットの計算もあらかじめ関数に落とし込んでおきます。ちょっとぶっ飛びますが、とりあえずコードを書きます。

#function2
def forward(X,T,params):
    W1,W2,B1,B2 = params
    H = sigmoid(np.dot(X,W1)+B1)
    Y = softmax(np.dot(H,W2)+B2)
    E = loss(Y, T)

    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)
    dE = [dW1, dW2, dB1, dB2]
    return [H, Y, E, dE]

まずは上から、「foward」という関数名でこの計算を定義して、そこには「X（入力）、T（正解ラベル）、prams（???）」の3つの変数を使いますよ、と記述します。このparamsですが、初期値をまとめたパラメータとして後の実際の計算のときに定義します。中身は「W1、W2、B1、B2」です（これ順番も大事）。

この関数の計算を始めるときには、paramsはそれぞれのパラメータが配列としてパックされた状態なので、使えるように展開してあげなければいけません。

W1,W2,B1,B2 = params

順番も大事ってのはここ、パックするパラメータの順番が間違ってると、ここで間違ったパラメータで展開されてしまいます。

このあとはよく知ってる計算が続きます。隠れ層（H）、出力層（Y）、誤差（E）まで計算して、微分式（dW1〜dB2）も合わせて計算しておきます。計算した微分式の結果は全部まとめてdEでパックしておきます。

dE = [dW1, dW2, dB1, dB2]

最終的に関数の計算はひとつしか出せないので、ひとつの配列としてほしい結果をすべてパックして出力します。

return [H, Y, E, dE]

こんな感じでひとつめの関数定義が終わり。

次にパラメータの更新も関数にしておきます、使う更新手法はモーメンタム法です。ちょっとややこしいかもしれないので、更新式を確認しておきます。

{\displaystyle 速度のパラメータ式：V_{new} = \gamma V_{old}-\mu \frac{\partial E}{\partial P} }

{\displaystyle パラメータ更新：P_{new} = P_{old}+V_{new} }

速度のパラメータ式と更新式、2つあるので、それぞれを関数で定義します。

まずは速度の式。

def velocity(learning_rate,momentum_term,Vs,dE):
    return [momentum_term*i-learning_rate*j for i, j in zip(Vs ,dE)]

2行ですが、とてもややこしい。まずは「velocity」という名前でこの関数を定義して、「learning_rate（学習率）、momentum_term（減衰係数）、Vs（更新したい速度のパラメータ：W1〜B2に対応）,dE（上で計算したやつ）」の5つのパラメータを使いますよ、と記述します。1行目でもうややこしい、ギリギリなんか分かる。

それで2行目、

{\displaystyle 速度のパラメータ式：V_{new} = \gamma V_{old}-\mu \frac{\partial E}{\partial P} }
この式がここに対応します。

momentum_term*i-learning_rate*j

iとjって何やねん、ってのをさらにその後ろのforで定義します。

for i, j in zip(Vs ,dE)

forって繰り返しのときに使うアレですよね、アレです、iとjはその後ろでzipされたVsとdEの中身にそれぞれ対応します。VsとdEはW1、W2、B1、B2、4つのパラメータに対応したものをまとめた変数なので、ひとつずつiとjで中身を取り出して計算する、という流れです。

次はパラメータの更新式、

def momentum(Vs,params):
    return [i+j for i, j in zip(params, Vs)]

ここも流れは同じです。paramsと上で計算したparamsを用いて「momentum」という名前の関数を定義します。

{\displaystyle パラメータ更新：P_{new} = P_{old}+V_{new} }
更新式はシンプルですね、iとjで計算するのがW1〜B2に対応したそれぞれparamsとVsの中身です。


なんということでしょう、この関数を使うことで、今まで8行使って書いていた計算を2行でまとめることができそうです。

〜ビフォア〜
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

〜アフタ〜
Vs = velocity(learning_rate,momentum_term,Vs,dE)
params = momentum(Vs,params)

関数の定義はこれで終了です。

初期値の設定

次に初期値の設定です。

#initialize
W1 = np.random.randn(2,3)
W2 = np.random.randn(3,2)
B1 = np.random.randn(1,3)
B2 = np.random.randn(1,2)
params_momentum = [W1,W2,B1,B2]

関数の中で使うparamsの設定です、テストで他の更新手法を試した経緯もあり、いちおう名前を「params_momentum」としています。ここでは名前にあんまり意味はないけど。

learning_rate = 0.008
momentum_term = 0.9
Vs = [np.zeros_like(i) for i in [W1,W2,B1,B2]]

速度のパラメータ設定です。初期値はすべてゼロにしているので、np.zerosを使って各パラメータの型に対応するVsをforで定義します。上の関数でやった形と同じですね。

save_momentum =[]

最後に損失関数を保存する配列を定義しておきます。

学習

ではイテレーションを回していきます。

#iteration
num_of_itr = 600
for i in range(num_of_itr):
    output_momentum = forward(X,T,params_momentum)
    Vs = velocity(learning_rate,momentum_term,Vs,output_momentum[3])
    params_momentum = momentum(Vs,params_momentum)

    #save loss
    save_momentum = np.append(save_momentum, output_momentum[2])

めっちゃシンプルですね、定義した関数を流れ通り使うだけです。各パラメータはリストでパックされているので、使うときには「output_momentum[3]」みたいにピックアップしてあげる必要があります。

グラフ化

ではグラフ化します。今回は隠れ層の結果がほしいのでHについて計算。

group1_hidden,group2_hidden = np.split(output_momentum[0],2)

まず、計算されたHは「output_momentum」でパックされているので「output_momentum[0]」で抜き取れます。上半分と下半分がそれぞれ分類したい2クラスに分かれているので、それをnp.splitします。

次にグラフ描画です。

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d') #3Dグラフのためのおまじない

3Dグラフにしたいので、おまじないが必要です。

その他、ラベルとか適当に設定します。

ax.set_xlabel('H1')
ax.set_ylabel('H2')
ax.set_zlabel('H3')
ax.scatter(group1_hidden[:,0], group1_hidden[:,1], group1_hidden[:,2])
ax.scatter(group2_hidden[:,0], group2_hidden[:,1], group2_hidden[:,2])
plt.show()

では結果を見てみましょう。
f:id:sutokun:20190318151012g:plain

ミスった、色をそろえてなかった、めんどくさいのでそのままいこう、すいません。青いドットセンターのクラスでオレンジのドットが外側のクラスに対応します。

これをみると最初はごちゃごちゃしていた2クラスのドットが3次元的に分離されていく様子が確認できます。へーすごい。

ちょっと見る角度を変えてみる、左がスタート、右が学習後。
f:id:sutokun:20190318151315p:plain:h200f:id:sutokun:20190318151337p:plain:h200

なんか学習後のグラフを見ていると、直線（平面）で分類できそうな気がします。つまり、隠れ層を入力層と見立てて、「入力層（隠れ層）→出力層」の関係で見てみると、これは3次元の線形分離問題そのもの。ごちゃごちゃのデータセットも隠れ層を入れることで分離され、非線形問題を線形分離問題に転換することができるんですね。

これが隠れ層の持つ大きな意味で、この隠れ層というブレイクスルーのおかげでディープラーニングは大きな発展を遂げることができたようです。

まとめ

今回は本題から少し外れて隠れ層のふるまいについて詳しく見てみました。なんで非線形問題には隠れ層が必要なのかというと、隠れ層が線形分類問題に転換してくれるんですね。ニューラルネットワークの良い理解に繋がりました。

全体コード
github.com

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

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

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

第二問

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

設問１．「原点を中心とした半径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

①非線形分類問題（２クラス分類）

第一問

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

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

ちょっと問題文がグダっていますね、どう書いていいのか分かりませんでした。要するに下図のようなデータセットを用意しろ、ということです。

f:id:sutokun:20190302134027j:plain:w450
中心のデータセットはこれまでと同じような作りでよさそうですね。問題はそれを取り巻く環状のデータセットです。問題と言っても円を書いてそこにガウシアンノイズをのせりゃいいだけなんですけどね。

ではコードで書いていきます。もちろん今回もワン・ホットで正解ラベルを書くのでライブラリのインポートにカラーマップはいりません。

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π)を使って書いてみましょう。

{\displaystyle 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

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

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

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

できあがった入力データセットの型はこんな感じ。

f:id:sutokun:20190302134121j:plain:w450

次に正解ラベルの作成です。

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を求めよ

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

f:id:sutokun:20190303005143j:plain:w450

隠れ層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の設定です。各層でのパラメータの型を確認してみます。
確認のために数式で書き下してみるとこんな感じ、ただしここでは行列を書きやすいように転置してるのでご注意。

f:id:sutokun:20190304021632j:plain:w450

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

#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:plainf:id:sutokun:20190304080214g:plain
何度かトライしてみると、「うまくいくとき」と「いかないとき」があります。体感では30％ぐらい失敗かな。

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

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

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

f:id:sutokun:20190304085119j:plain:w450
イメージを描いてみましたが、実際は多次元なのでこんなシンプルじゃないんですけどね。

この「ときどき損失量が収束しない問題」はディープラーニングを扱う以上どうしてもつきまといます。ただし、パラメータの更新を勾配降下法以外の更新方法を使うことで改善が見られます。次回は勾配降下法以外の更新方法でどう変わるのか見ていきたいと思います。


全体コード
github.com

次の問題
②非形分類問題（パラメータの更新） - 社畜エンジニア発掘戦線

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

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

④線形分類問題（多クラス分類）

第四問

前回はワン・ホット表現を用いて2クラス分類にトライしました。今回は更に多クラス（3クラス以上）の分類に取り組んでみたいと思います・

設問１．ガウシアンノイズを付加した入力x1、x2が(x1, x2) = (5, 5), (0, -5), (-5, 5)となるデータセットをワン・ホット表現を用いて作成せよ

クラス分類の入力が3つになりました。とは言え、前回の2クラスから3クラスになっただけなので、データグループの数を増やせばいいだけです。

まずはライブラリのインポート。

import numpy as np
import matplotlib.pyplot as plt

続いてデータセットを作成します。

#dataset
num_of_sam = 40
std_dv = 1.8

group1 = np.array([5,5])+np.random.randn(num_of_sam, 2)*std_dv
group2 = np.array([0,-5])+np.random.randn(num_of_sam, 2)*std_dv
group3 = np.array([-5,5])+np.random.randn(num_of_sam, 2)*std_dv
X = np.vstack((group1, group2, group3))

group3が増えていますが、今まで通りですね。では正解ラベルを書いていきます。

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

こちらも前回からgroup3が増えているだけですね。最後にXもTも結合するのを忘れないように注意です。

ではグラフ化してみましょう。

plt.scatter(group1[:,0],group1[:,1],marker='o',color='blue')
plt.scatter(group2[:,0],group2[:,1],marker='o',color='red')
plt.scatter(group3[:,0],group3[:,1],marker='o',color='green')
plt.show()

f:id:sutokun:20190227205352p:plain
確かに3クラスになっていますね、これでこの設問は終了です。

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

こちらも前回と流れは同じですが、データセットと出力の数が変わるので確認しておきます。
f:id:sutokun:20190227205424j:plain
線の数がだいぶ多くなってきました。

入力と出力の数は同じなので、基本的な骨格はこれまでと変わりません。なので、変更するのはパラメータW、Bの行列の型だけです。まずはソフトマックスとクロスエントロピーをPythonの関数に落としておきます。

#function
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
W = np.random.randn(2,3)
B = np.random.randn(1,3)

ニューロンの数が2→3となっているのでパラメータの型も変更しています。

では、準備が整ったので実際に計算してみます。

Y = softmax(np.dot(X, W)+B)
E = loss(Y, T)

ここは前回と同じですね。設問2は終了です。

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

パラメータの更新式ですが、行列の型が変わっているだけで式は同じです。

重み(W)の更新式：{\displaystyle W_{new} = W_{old} - \mu \Bigl( X \cdot (Y-T) \Bigr) }

バイアス(B)の更新式：{\displaystyle B_{new} = B_{old} - \mu \Bigl( 1 \cdot (Y-T) \Bigr) }

ではコードに落とします。学習率は0.003とします。

#initial setting
W = np.random.randn(2,3)
B = np.random.randn(1,3)

learning_rate = 0.003
E_save = []

#iteration
num_of_itr = 400
for i in range(num_of_itr):
    #forward propagation
    Y = softmax(np.dot(X, W)+B)
    E = loss(Y, T)
    E_save = np.append(E_save, E)

    #back propagation
    dW = X.T.dot(Y-T)
    dB = np.sum(Y-T, axis=0, keepdims=True)

    #update
    W = W - learning_rate*dW
    B = B - learning_rate*dB

ではグラフに描画します。グリッドを切ってそれぞれの出力を結合します。詳しくは前回の問題参照で。

#plot
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()]

Y_grid = softmax(np.dot(X_grid, W)+B)
Y_predict = np.around(Y_grid)

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

グリッドが切れたので実際にグラフ化します。

#plot_dataset
plt.scatter(group1[:,0],group1[:,1],marker='o',color='blue')
plt.scatter(group2[:,0],group2[:,1],marker='o',color='red')
plt.scatter(group3[:,0],group3[:,1],marker='o',color='green')

#plot_dataset
plt.scatter(green_group[:,0],green_group[:,1],marker='o',alpha=0.3,color='green')
plt.scatter(red_group[:,0],red_group[:,1],marker='o',alpha=0.3,color='red')
plt.scatter(blue_group[:,0],blue_group[:,1],marker='o',alpha=0.3,color='blue')
plt.show()

groupが3つになっていますが、書いている内容は同じですね。
f:id:sutokun:20190227205542p:plain
綺麗に分類できました、これで設問3は終了です。
[note] なぜか色の付かない欠損点ができるんですよね、出力ラベルが全部ゼロになってるのかな。

設問４．ガウシアンノイズを付加した入力x1、x2が(x1, x2) = (0, 0), (5, 5), (5, -5), (-5, -5), (-5, 5)となるデータセットを、勾配降下法で分類した結果と同一のグラフにプロットせよ

3クラス分類は2クラス分類からグループを増やしただけで実行することができたので、最後の設問として一気に5クラス分類にトライしてみたいと思います。

まずはデータセットの作成です。

import numpy as np
import matplotlib.pyplot as plt

#dataset
num_of_sam = 40
std_dv = 1.8

group1 = np.array([0,0])+np.random.randn(num_of_sam, 2)*std_dv
group2 = np.array([5,5])+np.random.randn(num_of_sam, 2)*std_dv
group3 = np.array([5,-5])+np.random.randn(num_of_sam, 2)*std_dv
group4 = np.array([-5,-5])+np.random.randn(num_of_sam, 2)*std_dv
group5 = np.array([-5,5])+np.random.randn(num_of_sam, 2)*std_dv
X = np.vstack((group1, group2, group3, group4, group5))

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

コピペコピペ＆コピペです。さすがに5グループ分も同じ記述が続くとブサイクに見えるので、この書き方は微妙です。もっとうまい記述はあると思いますが、とりあえず動きが見たいので今回コーディングの改良はスルーで。

グラフ化します。

plt.scatter(group1[:,0],group1[:,1], marker='o',color='yellow')
plt.scatter(group2[:,0],group2[:,1], marker='o',color='black')
plt.scatter(group3[:,0],group3[:,1], marker='o',color='red')
plt.scatter(group4[:,0],group4[:,1], marker='o',color='green')
plt.scatter(group5[:,0],group5[:,1], marker='o',color='blue')

plt.show()

f:id:sutokun:20190227210006p:plain
いい感じですね。

引き続きコピペを勧めます。活性化関数と損失関数は今まで通り「ソフトマックス」と「クロスエントロピー」です。

#function
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
W = np.random.randn(2,5)
B = np.random.randn(1,5)

learning_rate = 0.003
E_save = []

#iteration
num_of_itr = 400
for i in range(num_of_itr):
    #forward propagation
    Y = softmax(np.dot(X, W)+B)
    E = loss(Y, T)
    E_save = np.append(E_save, E)
    #back propagation
    dW = X.T.dot(Y-T)
    dB = np.sum(Y-T, axis=0, keepdims=True)
    #update
    W = W - learning_rate*dW
    B = B - learning_rate*dB

爆速で進んでいますが、ほぼ上記のコピペです。

さて、結果をグラフ化します。まずはグリッド切りから。

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

Y_grid = softmax(np.dot(X_grid, W)+B)
Y_predict = np.around(Y_grid)

out_connect = np.hstack((X_grid,Y_predict))
blue_group = out_connect[out_connect[:,2]==1]
green_group = out_connect[out_connect[:,3]==1]
red_group = out_connect[out_connect[:,4]==1]
black_group = out_connect[out_connect[:,5]==1]
yellow_group = out_connect[out_connect[:,6]==1]

実際にグラフ化します。

#plot_dataset
plt.scatter(group1[:,0],group1[:,1],marker='o',color='yellow')
plt.scatter(group2[:,0],group2[:,1],marker='o',color='black')
plt.scatter(group3[:,0],group3[:,1],marker='o',color='red')
plt.scatter(group4[:,0],group4[:,1],marker='o',color='green')
plt.scatter(group5[:,0],group5[:,1],marker='o',color='blue')

#plot_output
plt.scatter(blue_group[:,0],blue_group[:,1],marker='o',alpha=0.3,color='blue')
plt.scatter(red_group[:,0],red_group[:,1],marker='o',alpha=0.3,color='red')
plt.scatter(green_group[:,0],green_group[:,1],marker='o',alpha=0.3,color='green')
plt.scatter(black_group[:,0],black_group[:,1],marker='o',alpha=0.3,color='black')
plt.scatter(yellow_group[:,0],yellow_group[:,1],marker='o',alpha=0.3,color='yellow')

plt.show()

ほんと5グループも並ぶとコードはブサイクですね。まあ気にせず今回は力技でいきましょう。
f:id:sutokun:20190227210712p:plain
うまく分類できました。線形分類ですが5クラスちゃんと分けられるんですね。へーすごい。

まとめ

今回も学習が進む過程をGIFでアニメ化します、2本立てです。
f:id:sutokun:20190227210854g:plain　f:id:sutokun:20190227210828g:plain

さてさて、直線を使ってクラスを分類する「線形分類」を長々と扱ってきました。次回からは直線では分類できないようなデータセット（ドーナツ型とか）を分類できるような非線形問題にトライしていきたいと思います。これまではニューラルネットワークにとって比較的問簡単な問題だったので、パラメータの更新には勾配降下法しか使わなくてもそれなりに結果を出すことができました。しかし、そろそろ問題も複雑になってくるので、他の更新方法だとどのような違いが見られるのかも扱っていければと思います。

全体コード

  • 3クラス分類

github.com

  • 5クラス分類

github.com

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

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

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

③線形分類問題（２クラス分類：ワン・ホット エンコーダ）

第三問

今回はワン・ホット表現を用いた2クラスの分類問題にトライします。前回ちらっと書きましたが、2値（0と1）では2クラス以上の分類ができません（頭がいい人はできると思います）。そこで、3クラス以上の分類にも対応できるようにワンホット表現を導入して問題を解いていきます。

設問１．ガウシアンノイズを付加した入力x1、x2が(x1, x2) = (2.5, 2.5), (7.5, 7.5)となるデータセットをワン・ホット表現を用いて作成せよ

さて、そのワン・ホット表現ですが、もともとデジタル回路とかの分野で使われて表現のようです。

one-hot(ワン・ホット)は1つだけHigh(1)であり、他はLow(0)であるようなビット列のことを指す。
-by Wikipedia

2進数とはまたちょっと違うんですよね。例えば、5クラスをワンホットで表現すると、0=[1, 0, 0, 0, 0]、1=[0, 1, 0, 0, 0]、2=[0, 0, 1, 0, 0]、3=[0, 0, 0, 1, 0]、4=[0, 0, 0, 0, 1]、となります。プログラミングで数字は0から数えるので、5クラスだと最初は0で最後は4です。

今回の問題は2クラスなので入力x1の正解ラベルtは[1, 0]、x2は[0, 1]となります。では、ここまでをコードに落としていきます。

まずはライブラリのインポート、ワン・ホットではクラス別に直接色付けするので、今回カラーマップは必要ありません。

import numpy as np
import matplotlib.pyplot as plt

次にデータセットを作成します。

#dataset
num_of_sam = 40
std_dv = 1.8

group1 = np.array([2.5,2.5])+np.random.randn(num_of_sam, 2)*std_dv
group2 = np.array([7.5,7.5])+np.random.randn(num_of_sam, 2)*std_dv
X = np.vstack((group1, group2))

ここまでは今まで通りです。では正解ラベルを書いていきます。
[note]　上で書いたワン・ホットの順番が逆になっていますが、間違えただけです。順番が違うだけなので問題ないです。

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

イメージとしては[0, 1]をサンプル数だけnp.tileで並べていくって感じですかね。最後はTとして縦方向に結合します。

では、一度グラフ化してみます。2グループを連結したXとTは計算用です、グラフ化には連結前のそれぞれのデータセットを使います。ということは単にデータセットをグラフ化するだけなら正解ラベルTは必要なく、Xのデータセットに直接色付けをして描画します。

plt.scatter(group1[:,0],group1[:,1],marker='o',color='blue')
plt.scatter(group2[:,0],group2[:,1],marker='o',color='red')
plt.show()

f:id:sutokun:20190227104622p:plain
前回と同じようなグラフになりました。これでワン・ホット表現でのデータセットは完成です。

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

データセットができたところで、次はニューラルネットワークの設計です。ワン・ホット表現を用いるとクラスの数だけ正解ラベルと出力が必要になります。今回は2クラスなので出力は2つ、y1とy2を設定します。なので、設計するニューラルネットワークはこんな感じ。
f:id:sutokun:20190224032604j:plain

各パラメータをまとめて大文字の行列で記載しています。出力の数が変わっただけで骨格として大きな違いはありません。

さて新キャラのソフトマックス関数ですが、下記の式で表される関数です。

{\displaystyle y_{i} = \frac{e^{a_{i}}}{\sum_{j} e^{a_{j}}}\ }
iは計算が行われる出力ニューロンの順番を表しています。このソフトマックス関数は出力の「正解率（確率）」を計算してくれます。expがかかっていますが、分母のシグマは出力層すべてのニューロンの総和です。全体の総和を分母にして、各ニューロンを分子に置くことでそれぞれのニューロンからの出力を確率として表現してくれます。

このソフトマックス関数の特徴として、「出力の総和が1」となる性質があります。つまり、得られた出力がそのまま「正解率」として判断することができます。
f:id:sutokun:20190224035359j:plain

すごい便利ですね。ワン・ホットの正解ラベルは、「正解の順番の要素が1、それ以外は0」という表現でした。これは言い換えると、正解の確率は100%としてニューロンからの出力と比較することができます。こういう背景でワン・ホット表現とソフトマックス関数は相性が良いです。

では、このソフトマックス関数をPythonで関数としてコードに落としておきます。ついでにクロスエントロピーも合わせて書いておきます。

#function
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(2x2)とバイアスB(1x2)も定義しておきます。前回から出力の数が変わっているので設定する行列の型に注意です。

#initial setting
W = np.random.randn(2,2)
B = np.random.randn(1,2)

では計算してみます。

Y = softmax(np.dot(X, W)+B)
E = loss(Y, T)

関数を使って、そのままの記述ですね。ここまでで設問2は終了です。

設問３．勾配降下法でパラメータを400回更新し、その出力を求めよ

今回もいつ通り損失関数をパラメータで微分し、更新式を求めます。今回、パラメータは行列なので大文字でまとめて表記します。

重み(W)の微分：{\displaystyle \frac{\partial E }{\partial W}\ = \frac{\partial Y }{\partial W}\ \frac{\partial E }{\partial Y}\ }

バイアス(B)の微分：{\displaystyle \frac{\partial E }{\partial B}\ = \frac{\partial Y }{\partial B}\ \frac{\partial E }{\partial Y}\ }

で、今回は「ソフトマックス関数×クロスエントロピー」の計算なのでそれぞれの式は次のようになります。

重み(W)の微分：{\displaystyle \frac{\partial Y }{\partial W}\ \frac{\partial E }{\partial Y}\ = \frac{\partial }{\partial W}\ (X \cdot W + B) \cdot \Delta = X \cdot (Y-T) }

バイアス(B)の微分：{\displaystyle \frac{\partial Y }{\partial B}\ \frac{\partial E }{\partial Y}\ = \frac{\partial }{\partial B}\ (X \cdot W + B) \cdot \Delta = 1 \cdot (Y-T) }

この辺の計算はこちらを参照↓

以上より、勾配降下法を用いたパラメータの更新式は以下のようになります。

重み(W)の更新式：{\displaystyle W_{new} = W_{old} - \mu \Bigl( X \cdot (Y-T) \Bigr) }

バイアス(B)の更新式：{\displaystyle B_{new} = B_{old} - \mu \Bigl( 1 \cdot (Y-T) \Bigr) }

更新式が求められたので、これをコードに落とします。学習率は0.003とします。

#initial setting
W = np.random.randn(2,2)
B = np.random.randn(1,2)

learning_rate = 0.003
E_save = []

#iteration
num_of_itr = 400
for i in range(num_of_itr):
    #forward propagation
    Y = softmax(np.dot(X, W)+B)
    E = loss(Y, T)
    E_save = np.append(E_save, E)
    #back propagation
    dW = X.T.dot(Y-T)
    dB = np.sum(Y-T, axis=0, keepdims=True)
    #update
    W = W - learning_rate*dW
    B = B - learning_rate*dB

これで設問3は終了です。さすがに計算とそのコード書きはもう慣れてきましたね。

設問４．更新後のパラメータを用いた出力をデータセットと同じグラフ上にプロットせよ

更新式が求まったので、次はグラフへ出力します。グリッドを切ることろまではこれまでと同じです。

#plot_grid
grid_range = 10
resolution = 60
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()]

Y_grid = softmax(np.dot(X_grid, W)+B)
Y_predict = np.around(Y_grid)

ここからはいろいろやり方があると思うので、あくまで一例です（上級Python使いはどんなふうに書くのか気になります）。

まず、各グリッドごとの出力を横方向に結合します。この結合した行列から正解ラベル[1, 0]と[0, 1]に対応する行を抽出してblue_groupとred_groupに分けます、こんなイメージ。
f:id:sutokun:20190224074257j:plain

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

最後に、この分類したグループをグラフ化します。

#plot_dataset
plt.scatter(group1[:,0],group1[:,1],marker='o',color='blue')
plt.scatter(group2[:,0],group2[:,1],marker='o',color='red')

#plot_output
plt.scatter(red_group[:,0],red_group[:,1],marker='o',alpha=0.3,color='red')
plt.scatter(blue_group[:,0],blue_group[:,1],marker='o',alpha=0.3,color='blue')
plt.show()

f:id:sutokun:20190227104654p:plain
しっかり分類できていますね、これで設問4が終了です。

まとめ

今回も学習が進む過程をGIFでアニメ化します。
f:id:sutokun:20190227104713g:plain
ワン・ホットを使っても2値問題と同じように学習が進み、分類できることが確認できます。

今回の問題で2クラスをワン・ホットで分類できるようになりました。ワン・ホットが使えるようになったことで2クラス以上の分類も可能になります。次回の問題では多クラス分類にトライしてみようと思います。

全体コード
github.com

次の問題
④線形分類問題（多クラス分類） - 社畜エンジニア発掘戦線

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

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

活性化関数とクロスエントロピーの微分

おつかれさまです。

活性化関数や損失関数の微分についてメモしておきます。

〜活性化関数〜
ニューラルネットワークにおいて、線形変換をした後に適用する非線形関数もしくは恒等関数のこと。

よく使っている「シグモイド関数」が活性化関数ですね。

〜損失関数〜
ニューラルネットワークにおいて、出力(y)と正解(t)がどれだけ離れているのかを表す「損失量」を与える関数。

「二乗和誤差関数」、「クロスエントロピー」がこの損失関数にあたります。

今回は活性化関数に「シグモイド関数」、「ソフトマックス関数」、損失関数に「クロスエントロピー」を用いたときの微分の計算についてまとめてみたいと思います。
f:id:sutokun:20190218180627j:plain


まず、簡単なニューラルネットワークのモデルを下図のように設計します。今回の計算を分かりやすくするために入力(x)と出力(y)の間にaという中継変数を置きました。このaは単に入力とパラメータの線形結合を表す変数です。Eが損失関数、σは活性化関数を表しています。
f:id:sutokun:20190207094629j:plain
誤差逆伝搬法（バックプロパゲーション）では損失関数のパラメータ微分を行います。そして、計算が行えるように「連鎖律」というテクニックを用いて微分式を展開します。パラメータはwとbがありますが、とりあえずwで微分しましょう。中間変数aを入れたので連鎖律は三段構造になります。

{\displaystyle \frac{\partial E}{\partial w}\ = \frac{\partial a}{\partial w}\ \frac{\partial y}{\partial a}\ \frac{\partial E}{\partial y}\ }

活性化関数σにシグモイド関数を使う

では、この3つの微分をひとつずつ見ていきましょう。まず、ひとつ目の微分は線形結合wx+bを微分するだけなので簡単ですね。

{\displaystyle \frac{\partial a}{\partial w}\ = \frac{\partial }{\partial w}\ (xw+b) = x }

次の微分にはシグモイド関数を用いているので、その微分はσ(1-σ)で書けましたね。

{\displaystyle \frac{\partial y}{\partial a}\ = \frac{\partial }{\partial a}\ \sigma(a) = \sigma(a) (1-\sigma(a)) = y(1-y) }

最後に、クロスエントロピーの微分ですが、今回は入力xはたくさんの値を含んでいるわけではないので、シグマはいりません。

{\displaystyle E =\ -\sum_{i} \Bigl( t_{i} \log y_{i} +(1-t_{i}) \log (1-y_{i}) \Bigr) \rightarrow E = - ( t \log y +(1-t) \log (1-y) ) }

なので、その微分式は下記のようになります。

{\displaystyle \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial y}\ \Bigl( - ( t \log y +(1-t) \log (1-y) ) \Bigr) = -\frac{t}{y} + \frac{1 - t}{1 - y} }

ではもとの式に代入してみます。ただし、aのw微分の式はそのままにしておきます。

{\displaystyle \frac{\partial E}{\partial w}\ = \frac{\partial a}{\partial w}\ \frac{\partial y}{\partial a}\ \frac{\partial E}{\partial y}\ = \frac{\partial a}{\partial w}\ \cdot y(1-y)\cdot \Bigl(-\frac{t}{y} + \frac{1 - t}{1 - y} \Bigr) = \frac{\partial a}{\partial w}\ \cdot (y-t) }

この(y-t)をδと置きます。「なぜイキナリ新しい変数で置くねん」ってところはありますが、まずは置いてみます。

{\displaystyle \frac{\partial y}{\partial a}\ \frac{\partial E}{\partial y}\ =y(1-y)\cdot \Bigl(-\frac{t}{y} + \frac{1 - t}{1 - y} \Bigr) = (y-t) = \delta }

このδの意味ですが、「シグモイド関数×クロスエントロピー」の組み合わせで現れる微分結合式です。この結果を知っていれば、どれだけニューロンの数が増えても「中間変数aのパラメータw微分」だけ計算してあげればすぐに微分式が出せます。

つまり、

{\displaystyle \frac{\partial E}{\partial w}\ = \frac{\partial a}{\partial w}\ \delta = \frac{\partial }{\partial w}\ (wx+b) \cdot \delta = x \cdot (y-t) }

おースッキリしました。バイアスbで微分するときも同じです。

{\displaystyle \frac{\partial E}{\partial b}\ = \frac{\partial a}{\partial b}\ \delta = \frac{\partial }{\partial b}\ (wx+b) \cdot \delta = 1 \cdot (y-t) }

これはかなり計算が早くなりそうです。と言っても、出力層の計算なので隠れ層が何層もある場合はどれだけ効果が出るのか分かりませんが。

活性化関数σにソフトマックス関数を使う

分類問題によく出てくるソフトマックス関数は下記の式で表される関数です。

{\displaystyle y_{i} = \frac{e^{a_{i}}}{\sum_{j} e^{a_{j}}}\ }

ソフトマックス関数についてはここの問題でやや詳しめに書いています↓
③線形分類問題（２クラス分類：ワン・ホット エンコーダ） - 社畜エンジニア発掘戦線


こちらの関数も上記と同様に「クロスエントロピー」との結合微分式を求めたい、のですがめっちゃややこい。詳しい計算過程についてこちらのページに素晴らしく詳細に書かれていたので参照して下さい（あっ、コイツ逃げたぞ！とか言わないで…）。
qiita.com

とにかく（計算ができたとして）、微分結合式の結果はこのようになります。
※分類問題でソフトマックス関数を使う時は出力と正解ラベルが行列になることが多いので大文字で表記しています

{\displaystyle \frac{\partial Y}{\partial A}\ \frac{\partial E}{\partial Y}\ = }　（省略）　{\displaystyle = Y-T }

なんと、シグモイド関数と同じ結果となりました。つまり、ソフトマックス関数も損失関数にクロスエントロピーを使うならばΔ（大文字）はY-Tと表現できます。

{\displaystyle \frac{\partial E}{\partial W}\ = \frac{\partial A}{\partial W}\ \Delta = \frac{\partial }{\partial W}\ (W \cdot X + B) \cdot \Delta = X \cdot (Y-T) }

{\displaystyle \frac{\partial E}{\partial B}\ = \frac{\partial A}{\partial B}\ \Delta = \frac{\partial }{\partial B}\ (W \cdot X + B) \cdot \Delta = 1 \cdot (Y-T) }

これはシグモイド関数とソフトマックス関数、この2つの関数が深く関係していることに由来します。

まとめ

今回行った微分計算のイメージとしてはこんな感じ（分かりにくいけど）。
f:id:sutokun:20190218175934j:plain
活性化関数とクロスエントロピーの微分計算について少し詳しく（他のページに譲ったりして）見てみました。問題を解く上でこのδ（Δ）を使うこともあるかもしれないのでこの計算を参照してもらえばと思います。



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

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