pythonでニューラルネットワーク実装

# はじめに

pythonで3層のニューラルネットワークを実装し，XNORの識別をしてみました．
数式も載せたので，興味のある方は読んでみてください．
教科書として[『深層学習』](http://www.amazon.co.jp/深層学習-機械学習プロフェッショナルシリーズ-岡谷-貴之/dp/4061529021)を使いました．

本記事の構成

* はじめに
* ニューラルネットワーク
    * 重みの更新
    * 誤差逆伝播
* XNOR
* pythonでの実装
* 結果
* おわりに


# ニューラルネットワーク

ニューラルネットワークとは，人間の脳の神経回路を模したモデルです．このモデルを使うことで，画像認識や音声認識が可能となります．
今回実装したネットワークは入力層，中間層(1層)，出力層の3層構造です．

## 重みの更新

下の図で説明していきます．

<img width=360 alt="neu.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/8fb0bd49-0068-c66e-ebb8-c7c1e91cfbdf.png">

$l-1$ 層目の $i$ 番目のユニットから $l$ 層目の $j$ 番目のユニットへの重みを $w_{ji}^{(l)}$ とします．
また，$l-1$ 層目の $i$ 番目のユニットが保持している値を $u_{i}^{(l-1)}$ とします．
活性化関数としてシグモイド関数 $g(x)$ を用います．微分値が微分する前の値を使って求められるので便利です．

```math
g(x) = \cfrac{1}{1 + \exp(-x)} \\
g'(x) = g(x)(1 - g(x))
```

2乗誤差 $E_n$ を目的関数とし，下の式で重みを更新することで最小化を図ります．$E_n$ は1サンプルによって生じる誤差で，この $E_n$ を使って重みの更新をする方法を確率的勾配降下法と呼びます．
確率的勾配降下法は，目的関数が更新されるごとに変化するため，局所解にトラップされにくくなるというメリットがあります．
$\epsilon$ は学習率と呼ばれ，学習の速度を決めるパラメータです．
偏微分の項をどのように計算するかがポイントとなります．

```math
{{w}_{ji}^{(l)}}_{new} = {w_{ji}^{(l)}}_{old} - \epsilon\cfrac{\partial E_n}{\partial {w_{ji}^{(l)}}_{old}}
```

## 誤差逆伝播

では，$\cfrac{\partial E_n}{\partial w_{ji}^{(l)}}$ の求め方を説明していきます．出力層と中間層とで求め方が若干異なるので，分けて説明していきます．

**< 出力層 >**
教師を $t_j$ で表します．2乗誤差 $E_n$ を出力層の重み $w_{ji}^{(l)}$ で偏微分すると，

```math
\begin{align}
\cfrac{\partial E_n}{\partial w_{ji}^{(l)}} &= \cfrac{\partial}{\partial w_{ji}^{(l)}}\cfrac{1}{2}(t_j - g(u_{j}^{(l)}))^{2} \\
&= (t_j - g(u_{j}^{(l)}))\cdot\cfrac{\partial}{\partial w_{ji}^{(l)}}(t_j - g(u_{j}^{(l)})) \\
&= (t_j - g(u_{j}^{(l)}))\cdot\cfrac{\partial}{\partial u_{j}^{(l)}}(t_j - g(u_{j}^{(l)}))\cdot\cfrac{\partial u_{j}^{(l)}}{\partial w_{ji}^{(l)}} \\
\\
&= (g(u_{j}^{(l)}) - t_j)\cdot g'(u_{j}^{(l)})\cdot g(u_{i}^{(l-1)}) \\
\\
&= (g(u_{j}^{(l)}) - t_j)g(u_{j}^{(l)})(1 - g(u_{j}^{(l)}))g(u_{i}^{(l-1)})
\end{align}
```

**< 中間層 >**
2乗誤差 $E_n$ を中間層の重み $w_{ji}^{(l)}$ で偏微分すると，

```math
\begin{align}
\cfrac{\partial E_n}{\partial w_{ji}^{(l)}} &= \cfrac{\partial E_n}{\partial u_{j}^{(l)}}\cfrac{\partial u_{j}^{(l)}}{\partial w_{ji}^{(l)}} \\
\\
&= \delta_{j}^{(l)}\cfrac{\partial u_{j}^{(l)}}{\partial w_{ji}^{(l)}}
\end{align}
```

右辺第1項は，

```math
\begin{align}
\delta_{j}^{(l)} &= \cfrac{\partial E_n}{\partial u_{j}^{(l)}} \\
&= \sum_{k}\cfrac{\partial E_n}{\partial u_{k}^{(l+1)}}\cfrac{\partial u_{k}^{(l+1)}}{\partial u_{j}^{(l)}} \\
\\
&= \sum_{k}\delta_{k}^{(l+1)}w_{kj}^{(l+1)}g'(u_{j}^{(l)}) \\
\\
&= \Bigl(\sum_{k}\delta_{k}^{(l+1)}w_{kj}^{(l+1)}\Bigr)g(u_{j}^{(l)})(1 - g(u_{j}^{(l)}))
\end{align}
```

右辺第2項は，

```math
\begin{align}
\cfrac{\partial u_{j}^{(l)}}{\partial w_{ji}^{(l)}} &= g(u_{i}^{(l-1)})
\end{align}
```

以上より，

```math
\begin{align}
\cfrac{\partial E_n}{\partial w_{ji}^{(l)}} &= \delta_{j}^{(l)}\cfrac{\partial u_{j}^{(l)}}{\partial w_{ji}^{(l)}} \\
\\
&= \Bigl(\sum_{k}\delta_{k}^{(l+1)}w_{kj}^{(l+1)}\Bigr)g(u_{j}^{(l)})(1 - g(u_{j}^{(l)}))g(u_{i}^{(l-1)})
\end{align}
```

中間層 $l$ の重みを更新するときは，その次の層 $l+1$ の $\delta_{k}^{(l+1)}$ が必要となります．
出力層での $\delta$ は，出力値と教師の差として求めることができ，これを出力層から入力層の方向に順番に伝播させることで，中間層の重みを更新します．これが誤差逆伝播と呼ばれる所以です．


# XNOR

XNORは入力値が等しいとき $1$ を出力し，入力値が異なるとき $0$ を出力します．これは線形識別することができません．

|$x_1$|$x_2$|$t$|
|-----|-----|---|
|0    |0    |1  |
|0    |1    |0  |
|1    |0    |0  |
|1    |1    |1  |

<img width=360 alt="xnor.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/afa16cf9-6af1-63ef-13fe-3306fd54f5b4.png">


# pythonでの実装

XNORを識別するニューラルネットワークを実装しました．
中間層の数は1つ，中間層のユニットの数は2つにしています．学習率 $\epsilon = 0.1$，モメンタムの係数 $\mu = 0.9$ としています．
(モメンタムは収束性能を向上させる方法の一つで，重みの修正量に，前回の重みの修正量に係数をかけたものを加算します．)

```neuralnetwork.py
import numpy
import math
import random
from matplotlib import pyplot

class Neural:

	# constructor
	def __init__(self, n_input, n_hidden, n_output):
		self.hidden_weight = numpy.random.random_sample((n_hidden, n_input + 1))
		self.output_weight = numpy.random.random_sample((n_output, n_hidden + 1))
		self.hidden_momentum = numpy.zeros((n_hidden, n_input + 1))
		self.output_momentum = numpy.zeros((n_output, n_hidden + 1))


# public method
	def train(self, X, T, epsilon, mu, epoch):
		self.error = numpy.zeros(epoch)
		N = X.shape[0]
		for epo in range(epoch):
			for i in range(N):
				x = X[i, :]
				t = T[i, :]

				self.__update_weight(x, t, epsilon, mu)

			self.error[epo] = self.__calc_error(X, T)


	def predict(self, X):
		N = X.shape[0]
		C = numpy.zeros(N).astype('int')
		Y = numpy.zeros((N, X.shape[1]))
		for i in range(N):
			x = X[i, :]
			z, y = self.__forward(x)

			Y[i] = y
			C[i] = y.argmax()

		return (C, Y)


	def error_graph(self):
		pyplot.ylim(0.0, 2.0)
		pyplot.plot(numpy.arange(0, self.error.shape[0]), self.error)
		pyplot.show()


# private method
	def __sigmoid(self, arr):
		return numpy.vectorize(lambda x: 1.0 / (1.0 + math.exp(-x)))(arr)


	def __forward(self, x):
		# z: output in hidden layer, y: output in output layer
		z = self.__sigmoid(self.hidden_weight.dot(numpy.r_[numpy.array([1]), x]))
		y = self.__sigmoid(self.output_weight.dot(numpy.r_[numpy.array([1]), z]))

		return (z, y)

	def __update_weight(self, x, t, epsilon, mu):
		z, y = self.__forward(x)

		# update output_weight
		output_delta = (y - t) * y * (1.0 - y)
		_output_weight = self.output_weight
		self.output_weight -= epsilon * output_delta.reshape((-1, 1)) * numpy.r_[numpy.array([1]), z] - mu * self.output_momentum
		self.output_momentum = self.output_weight - _output_weight

		# update hidden_weight
		hidden_delta = (self.output_weight[:, 1:].T.dot(output_delta)) * z * (1.0 - z)
		_hidden_weight = self.hidden_weight
		self.hidden_weight -= epsilon * hidden_delta.reshape((-1, 1)) * numpy.r_[numpy.array([1]), x]
		self.hidden_momentum = self.hidden_weight - _hidden_weight


	def __calc_error(self, X, T):
		N = X.shape[0]
		err = 0.0
		for i in range(N):
			x = X[i, :]
			t = T[i, :]

			z, y = self.__forward(x)
			err += (y - t).dot((y - t).reshape((-1, 1))) / 2.0

		return err
```


```main.py
from neuralnetwork import *

if __name__ == '__main__':

	X = numpy.array([[0, 0], [0, 1], [1, 0], [1, 1]])
	T = numpy.array([[1, 0], [0, 1], [0, 1], [1, 0]])
	N = X.shape[0] # number of data

	input_size = X.shape[1]
	hidden_size = 2
	output_size = 2
	epsilon = 0.1
	mu = 0.9
	epoch = 10000

	nn = Neural(input_size, hidden_size, output_size)
	nn.train(X, T, epsilon, mu, epoch)
	nn.error_graph()

	C, Y = nn.predict(X)

	for i in range(N):
		x = X[i, :]
		y = Y[i, :]
		c = C[i]

		print x
		print y
		print c
		print ""
```


# 結果

正しく識別できています．

|$x_1$|$x_2$|$t$      |$y$                        |
|:---:|:---:|:-------:|:-------------------------:|
|0    |0    |[ 1, 0 ] |[ 0.92598739, 0.07297757 ] |
|0    |1    |[ 0, 1 ] |[ 0.06824915, 0.93312514 ] |
|1    |0    |[ 0, 1 ] |[ 0.06828438, 0.93309010 ] |
|1    |1    |[ 1, 0 ] |[ 0.92610205, 0.07220633 ] |

誤差が減少していく様子をグラフにしました．横軸がエポック数，縦軸が誤差です．

<img width=840 alt="error.png" src="https://qiita-image-store.s3.amazonaws.com/0/82527/865e697e-2f9b-69ce-e579-e5b9e1e2e723.png">


# おわりに

3層のニューラルネットワークを実装し，XNORを識別することができました．
学習のトリック(重みの初期値の決定方法，学習率の決定方法，AdaGrad，モメンタム)など勉強すべきことはまだまだありますが，一旦これで終わりにします．