社畜エンジニア発掘戦線

駆けだしAIエンジニア

①非線形回帰問題（データフィッティング）

第一問


直線をデータにフィットさせる線形問題はいくつか取り組んできました。ここからは直線ではフィットできないデータにカーブをフィットさせる非線形問題にチャレンジしていきたいと思います。

今回もPYTHONISTA3でコードを書きながら解いていってみましょう。

設問１．ガウシアンノイズを付加したy = sin(x)に準ずるデータセット（インプット：x、正解ラベル：t）を作成せよ

さて、今回から非線形問題ということで、データセットする関数が直線ではなくsinカーブとなりました。とは言え、基本的に解決フローは今までと同じです。そしてデータセットは直線：y=oxのところがy=sin(x)に置きかわるだけです。範囲はsinのだいたい1周期にあたる0から6まででトライしましょう。サンプル数は40にします。ノイズの広がりが大きいとフィッティングが乱れやすいので、標準偏差は0.1で極力狭くしておきます。

import numpy as np
import matplotlib.pyplot as plt

num_of_sam = 40
x = np.random.uniform(0, 6, num_of_sam)
noise = np.random.normal(0, 0.1, num_of_sam)
t = np.sin(x)+noise

グラフにプロットしてみます。

plt.plot(x, t, 'o')
plt.show()

f:id:sutokun:20190204062421p:plain
オッケーですね、設問１は終了です。

設問２．隠れ層が20のニューロンで設計されるニューラルネットワークのパラメータ（w、b）の行列型を求めよ

ココから線形問題よりも難易度がぐっと上がります。ですが、解決フローに則って進めていくことに変わりはないので、ひとつずつトライしていきましょう。

まずは、ニューラルネットワークの設計です。非線形と名のつく問題には隠れ層が必須です。入力と出力の間に位置する隠れ層が20のニューロンで構築されるので、次のような感じになります。ニューロンとは矢印を受けたり出したりする玉のことです。
f:id:sutokun:20190127151059j:plain
いきなりめっちゃパラメータが増えましたね。入力（x）から入るデータセットは隠れ層（h）を通り、出力（y）として出ていきます。隠れ層で実行される計算は今まで通りxw+bですが、ここに非線形関数σを寄与させなければいけません。つまり、隠れ層の出力はh=σ(wx+b)となります。

なぜ非線形関数（曲線）を寄与させなければいけないのかというと、早い話が「直線をいくら足してもいいフィッティングにはならない」からです。wx+bのままでは直線ですが、σに入れてやることで曲線に変換できます。そんな感じの理由です。

ここに用いる関数も場面によって使い分ける必要がありますが、スタンダードなのはシグモイド関数と呼ばれるものです。

この関数のグラフはこんな形をしています。
f:id:sutokun:20190127151027j:plain

このシグモイド関数を２０個（ニューロン数）足し合わせていい感じのフィッティングカーブを作ろう！というのが大まかな流れです。ちなみに出力yには特になんの関数も必要なく線形結合（ただの足し合わせ）で大丈夫です。

では、シグモイド関数をPythonの関数に落としておきましょう。

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

これで「sigmoid(oo)」とするだけでシグモイド関数の計算ができます。

ひと通り準備も終えたところで、パラメータについて考えてみます。隠れ層のニューロン数は20なので重み（w）とバイアス（b）はそれぞれ20個必要です。次に、隠れ層からの出力20個が、ニューラルネットワーク全体の出力（y）の入力となります。それぞれに重み（w）が掛けられて出力（y）へ向かいますが、出力に作用するバイアスはひとつだけです。
f:id:sutokun:20190127155152j:plain

整理すると、最初のwとbは(20, 1)のデータとなり、その次はwが(1, 20)、bは(1, 1)、つまりスカラー値となります。wは前と後ろで型が反転していることに注意です。出力は線形結合なので反転させて内積計算にしなければいけません。

コードに落としましょう。W1とW2が前側のパラメータなので(20, 1)となり、W2とB2が後ろ側のパラメータなのでW2は(1, 20)、B2は(1, 1)となっています。いちいち全部の初期値を設定していると面倒なので、ランダムで値が入るようにしています。

dimention = 20
W1 = np.random.randn(dimention, 1)
B1 = np.random.randn(dimention, 1)

W2 = np.random.randn(1, dimention)
B2 = np.random.randn(1, 1)

これで設問２が終了です。

設問３．ニューラルネットワークの出力（y）を計算せよ

パラメータの初期値も設定できたので、実際に計算を進めていきたいと思います。まずは隠れ層の出力です。

H = sigmoid(x*W1+B1)

型を確認しておきましょう。入力（x）は40個のデータなので40列のデータがニューロン数の20行に収まっています。

print(H.shape) #(20,40)

次は、計算したHを用いて出力yを求めます。線形結合なので内積にしなければいけません。dotは内積の意味です。

y = np.dot(W2, H)+B2

こちらも念のため、型が正しいか確認しておきましょう。入力（x）と一致するはずです。が、この計算だとyは行列として出力されますが、もとのxは配列です。そのためshapeで出力される型は異なります。これでも数値計算は進みますが、最後グラフ化するときに型をそろえる必要があります。

print(x.shape) #(40, )
print(y.shape) #(1, 40)
設問４．二乗和誤差関数から損失量を求め、各パラメータによる損失関数の微分式を導出せよ

次は損失関数です、今回も二乗和誤差関数を用いましょう。

E = np.sum((t-y)**2)

さて、ここからは計算です。損失関数を４つのパラメータ（W1、B1、W2、B2）で微分します。これは逆伝搬計算にあたるのでW2、B2から計算していきましょう。その前に、Eをyで微分した結果はよく使うので先にこれだけ計算して変数δで置いておきます。

{\displaystyle \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial y}\ \sum_{i} (t_{i} - y_{i})^2 = 2 \sum_{i} (y_{i}-t_{i}) = \delta_{i} }

まずはW2、B2から計算です、これまで同様「連鎖率」を使います。

{\displaystyle \frac{\partial E}{\partial W_2}\ = \frac{\partial y}{\partial W_2}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial W_2}\ \Bigl(\sum_{i} (H_{i}W_2)+B_2\Bigr) \cdot \delta_{i} = 2\sum_{i} H_{i} (y_{i}-t_{i}) }
{\displaystyle \frac{\partial E}{\partial B_2}\ = \frac{\partial y}{\partial B_2}\ \frac{\partial E}{\partial y}\ = \frac{\partial }{\partial B_2}\ \Bigl(\sum_{i} (H_{i}W_2)+B_2\Bigr) \cdot \delta_{i} = 2\sum_{i} (y_{i}-t_{i}) }

次にW1、B2の計算です。しかし、その前にシグモイド関数の微分を確認しておきましょう。

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

ややこしい中身の計算は紙とペンでトライしてもらうとして、結果だけ書くと非常にシンプルな式にまとまります。

このシグモイド関数の微分結果を用いて1層目のパラメータ微分を進めます。こちらも同様プロセスですが、2層分戻らなければいけないのでダブルで「連鎖律」を使います。

{\displaystyle \frac{\partial E}{\partial W_1}\ = \frac{\partial H}{\partial W_1}\ \frac{\partial y}{\partial H}\ \frac{\partial E}{\partial y}\ = \frac{\partial H}{\partial W_1}\ \frac{\partial y}{\partial H}\ \cdot \delta_{i} }

ボリュームがあるのでひとつずつ計算します。

{\displaystyle \frac{\partial H}{\partial W_1}\ = \frac{\partial }{\partial W_1}\ \sigma (x_{i}W_1+B_1) = \sum_{i} x_{i}H_{i}(1-H_{i}) }
{\displaystyle \frac{\partial y}{\partial H}\ = \frac{\partial }{\partial H}\ (H_{i}W_2+B_2) = W_2 }

以上をまとめると、

{\displaystyle \frac{\partial E}{\partial W_1}\ = W_2 \sum_{i} x_{i}H_{i}(1-H_{i}) \cdot \delta_{i} =2 W_2 \sum_{i} x_{i}H_{i}(1-H_{i})(y_{i}-t_{i}) }

こりゃ大変ですね、B1については結果だけ。

{\displaystyle \frac{\partial E}{\partial B_1}\ = W_2 \sum_{i} H_{i}(1-H_{i}) \cdot \delta_{i} =2 W_2 \sum_{i} H_{i}(1-H_{i})(y_{i}-t_{i}) }

なんとか微分式の導出までできたので、コードに落としておきます。

dW2 = 2*np.sum(H*(y-t),axis=1)
dB2 = 2*np.sum(y-t)

dW1 = 2*W2*np.sum(x*H*(1-H)*(y-t),axis=1)
dB1 = 2*W2*np.sum(H*(1-H)*(y-t))

dW2、dW1、dB1はHの項が絡むのでニューロンの数だけ個別にパラメータの更新が必要です。そのため、axis=1を加えて行要素のみにsumが効くように調整しています。…なんですが、バイアスの変化はとてもシビアに出力へ影響します。バイアスを個別に更新すると損失関数がうまく収束に向かわない傾向があるのでdB1からはaxis=1を除いて、全ての要素のsumで更新されるようにしています。

長くなりましたが、これで設問４が終了です。

設問５．降下法を用いてパラメータの更新式を導出せよ

微分式が求められたのでココはもう簡単ですね。Pはパラメータ（w、b）の意味です。

降下法：{\displaystyle P_{new} = P_{old} - \mu \frac{\partial E}{\partial P}\ }

この式に先の微分式を代入するだけ。

{\displaystyle W_{2new} = W_{2old} - \mu \frac{\partial E}{\partial W_{2}}\ = W_{2old} - \mu \Bigl(2 W_2 \sum_{i} x_{i}H_{i}(1-H_{i})(y_{i}-t_{i})\Bigr) }
{\displaystyle B_{2new} = B_{2old} - \mu \frac{\partial E}{\partial B_{2}}\ = B_{2old} - \mu \Bigl(2 W_2 \sum_{i} H_{i}(1-H_{i})(y_{i}-t_{i})\Bigr) }

{\displaystyle W_{1new} = W_{1old} - \mu \frac{\partial E}{\partial W_{1}}\ = W_{1old} - \mu \Bigl( 2\sum_{i} H_{i} (y_{i}-t_{i})\Bigr) }
{\displaystyle B_{1new} = B_{1old} - \mu \frac{\partial E}{\partial B_{1}}\ = B_{1old} - \mu \Bigl( 2\sum_{i} (y_{i}-t_{i})\Bigr) }

ではコードに落とします。学習率（learning_rate）はこの後定義するので、とりあえずコードを書いてしまいましょう。

W1 = W1-learning_rate*dW1.T
W2 = W2-learning_rate*dW2
B1 = B1-learning_rate*dB1.T
B2 = B2-learning_rate*dB2

W1とB1は計算される出力は同じなのですが、型が反転します。そのためTで転置させています。（この辺は本質ではないのでコーディングが上手な方はうまくアレンジして下さい）

設問５が終了です。

設問６．学習率を0.002として、イテレーションを3000回繰り返した後の出力を求めよ

ようやくゴールが見えてきました。イテレーションの数がハンパなく増えましたが、非線形の収束にはこれくらい必要なようです。ではいつもも通りforを使って繰り返しのコードを書いてみましょう。どうせこの後にコスト関数のグラフ化をするので、先に損失量を保存するコードも入れておきます。

learning_rate = 0.002
E_save = []

num_of_itr = 3000
for i in range(num_of_itr):
    #forward propagation
    H = sigmoid(x*W1+B1)
    y = np.dot(W2,H)+B2
    E = np.sum((t-y)**2)
    E_save = np.append(E_save,E)

    #differential
    dW2 = 2*np.sum(H*(y-t),axis=1)
    dB2 = 2*np.sum(y-t)
    dW1 = 2*W2*np.sum(x*H*(1-H)*(y-t),axis=1)
    dB1 = 2*W2*np.sum(H*(1-H)*(y-t))

    #back propagation
    W1 = W1-learning_rate*dW1.T
    W2 = W2-learning_rate*dW2
    B1 = B1-learning_rate*dB1.T
    B2 = B2-learning_rate*dB2

では、ご待望のグラフ化です。描画するのはカーブなのでx軸は200点と細かめに刻んでおきます。最終的に出力されたW1、B1、W2、B2を使って描画用の入力（X_line）を順伝搬計算します。が、このままではグラフが描けません。上の方でちょろっと書いたんですが、この順伝搬計算の出力は行列で出力されますが、入力は配列です。そのためそれらの型をそろえないとグラフが描けません。出力はravelを使って配列にしてしまいます。（ここも本質ではないので、コーディングが上手な方はアレンジして下さい）

X_line = np.linspace(0, 6, 200)
H_line = sigmoid(X_line*W1+B1)
Y_line = np.ravel(np.dot(W2,H_line)+B2)

plt.plot(x, t, 'o')
plt.plot(X_line, Y_line, 'o')
plt.show()

f:id:sutokun:20190204062717p:plain
尻尾のあたりが微妙ですが頑張ってフィッテイングしてくれました感じは伝わります。

設問７．損失量のイテレーション回数による推移をグラフ化せよ

いつも通り、損失関数で求められる損失量を確認してみましょう。既に損失量を保存するコードは書いているので、それをプロットするだけです。

plt.plot(E_save)
plt.show()

f:id:sutokun:20190204062803p:plain
大きな問題なく、収束に向かってくれています。尻尾のあたりがうまくフィットしないのに損失量の推移は落ち着いてしまってこれ以上フィッティングが進みません。んん、これはどうしたもんかな〜という感じです。

まとめ

非線形問題の第一問が完了しました。同じくGIFでアニメ化します。イテレーションは3000回でそんなに画像出力はできないので、100回ごとに画像出力させて全部で30枚をアニメ化しています。
f:id:sutokun:20190204063008g:plain

大きな問題解決フローは今までと変わりませんが、なにせパラメータの数が一気に増えたので頭が爆発しそうになります。ニューラルネットワークを使った非線形フィッティングはあくまでシグモイド関数の足し合わせで表現されるので、完ぺきにきれいな曲線になることはありません。じゃあシグモイド関数以外の関数を使ったらどうなるんだ、とか面白そうです。また問題にして考えてみたいと思います。隠れ層の出力をグラフ化してシグモイド関数たちがどう動くのか、とかも面白そう。

全体コード
github.com

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

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

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