多層パーセプトロンの実装(Python)

今回は自分の勉強を兼ねて、順伝播型ニューラルネットワークであるパーセプトロンの基礎知識と、numpyを使った実装についてまとめました。

パーセプトロン(順伝播型NN)の構造

教師あり学習のタスクにおいては、どのような入力を与えるかが学習の精度に大きく関わってきます。

データの中には人の手によって選別できるような特徴量(その現象にとって本質的な情報)もありますが、より高い精度を追求する場合、人の手による選別では不十分になってきます。

そこで、特徴量の抽出も機械にやらせてしまおうということで生まれたのがニューラルネットです。

ニューラルネットは入力層〜中間層でタスクに適したデータの表現を獲得し、出力層で回帰等のタスクを実行するという流れで推論が行われます。

説明の前にまず記号を定義しておきます。

$$\begin{align}z^{(l)}_i&:第 l層の j番目のユニットの出力\\
w^{(l)}_{ji}&:第 l層の重みベクトルの ji成分\\
b^{(l)}_j &:第 l層に固有のバイアス項の第j成分\\
u_j^{(l)}=\sum_{i} w^{(l)}_{ji}z^{(l-1)}_i&: 第 l 層の j 番目のユニットに対する入力\end{align}$$

入力層

パーセプトロンへの入力ベクトルの各成分を出力するだけのユニットです。

入力ベクトルの次元数と等しいユニットを持ち、活性化関数は持ちません。(強いていうなら恒等写像が活性化関数)

式に直すと、

$$z_i^{(1)}=x_i, z^{(1)} = x\\$$

(ただし、入力を$x=(x_1,x_2,…,x_n)^T$とする)

となります。

中間層

入力データからそのタスクに適したデータの表現を獲得する役割をもっています。

第$l$層の活性化関数を$f^{(l)}$とすると、第$l$層の出力は、$z^{(l)}_j=f^{(l)}(b^{(l)}_j+u_j^{(l)})$です。

$u^{(l)}=(u_j^{(l)}),z^{(l)}=(z^{(l)}_j),b^{(l)}=(b_j^{(l)}),W=\begin{pmatrix}w^{(l)}_{11}&w^{(l)}_{12}&\dots\\w^{(l)}_{21}&w^{(l)}_{22}&\dots\\.& &.\\.&&&.\\\end{pmatrix}$とおいて中間層の演算をベクトル表記すると、以下のようになります。

$$u^{(l)}=W^{(l)}z^{(l)}, z^{(l)}=f^{(l)}(u^{(l)}+,b^{(l)})$$

ここで、$f^{(l)}$はベクトルの各成分に作用するものとします。

今回は活性化関数$f^{(l)}$としてシグモイド関数$\frac{1}{1+e^{-x}}$を使用しました。

出力層

出力層は、回帰や分類などのタスクをこなす部分になります。

そのためバイアス項は使用しません。

今、パーセプトロンが$L-1$層からなるとします。この時の出力は以下のようになります。

$$\hat{y}=z^{(L)}=f^{(L)}(u^{(L)}),\quad u^{(L)}=W^{(L)}h(ただしh = z^{(L-1)})$$

以上をまとめると、順伝播型ニューラルネットとは、入力$x$の表現$h=z^{(L-1)}$の回帰分析を通じて$x$と$y$の関係を推定する仕組みである、ということになります。

最適化手法:勾配降下法(gradient descent)

一般的に機械学習のタスクは、入力をもとに損失関数の値を最小化することで尤もらしい推論を行います。

線形回帰の場合はMSE(mean squared estimator)、分類問題の場合はクロスエントロピー(分布関数の対数尤度関数)を使用することが多いです。

いずれの場合も連続関数の最小問題に帰着されるため、数値計算によってその解を求めることになります。

計算コストの関係上、一階微分の情報をもとにして近似的に計算をしていくことになります。
(Newton法のように二階微分の項が出てくるとヘッセ行列の逆行列計算が必要になり計算コストが嵩む。O($N^3$)なので現実的でない。)

そこでよく用いられるのが勾配降下法という手法です。

アイデアとしては、勾配が0になるまで損失関数の勾配方向に移動を続けることで極小値を見つけよう、というイメージです。

損失関数を$E(w)$、求めたい最適解を$w^*=argmin_{w}E(w)$とします。

なるべく低い位置(目的関数値の小さいところ)に進みたいので、重みベクトルを目的関数の勾配を逆方向に動かしていけば良いです。時間パラメータ$t$に対して離散化した定式化は以下のようになります。

$$w^{(t+1)}=w^{(t)}+\Delta w^{(t)},\Delta w^{(t)} = -\eta \nabla E(w^{(t)})$$

ここで$\eta$は学習率と呼ばれるハイパーパラーメーターで、これを大きくすると学習が早く進みますが目的関数値が収束しない可能性があります。小さくしすぎると学習に時間がかかり過ぎてしまうため、適切な値を使用することが重要です。

勾配降下法は近似解法ですので、初期値の選び方によって暫定解が変化したり、解がうまく収束しなかったりと問題点もあります。

それらを改善するために改良されたアルゴリズムもたくさんありますが、今回は実装していませんので割愛させていただきます。

次のセクションでは重みベクトルの初期値の選び方について紹介します。

重みベクトルの初期化方法

LeCunの初期化

分散$\frac{1}{d_{l-1}}$の一様分布またはガウス分布からサンプリングする方法。($d_{l-1}$は第$l-1$層のユニット数)

$w_{ji}^{(l)}$ ~ $\mathscr{U}(w;-\sqrt{\frac{3}{d_{l-1}}} , \sqrt{\frac{3}{d_{l-1}}})$ or $ \mathcal{N}(w;0,\frac{1}{\sqrt{d_{l-1}}})$

このようにサンプリングすることで、前層との結合を多く持つユニットに対しては重みの初期値を小さく抑えることができます。

一方でCNNなど疎なモデルの場合には第$l$層のユニットと結合している前層のユニットの総数をfan-inとして新たに定義し、$d_{l-1}$の代わりに使用します。

Glorotの初期化

線形ユニットのみを持つNNの解析から生まれた手法。左右対称な形をした活性化関数に対して有用です。なぜ、どのように有用なのかはこれから勉強します。

サンプリング方法は以下の通り。

$w_{ji}^{(l)}$ ~ $\mathscr{U}(w;-\sqrt{\frac{6}{d_{l-1}+d_{l}}} , \sqrt{\frac{6}{d_{l-1}+d_{l}}})$ or $ \mathcal{N}(w;0,\sqrt{\frac{2}{d_{l-1}+d_{l}}})$

一般のネットワークの場合、$d_{l-1}$をfan-inに、$d_{l}$をfan-out(次の層のユニットのいくつと結合しているか)に置き換えて使用します。

Heの初期化

ReLUを用いたNNの解析の過程で生まれた手法です。サンプリング手法は以下の通り。

$w_{ji}^{(l)}$ ~ $\mathscr{U}(w;-\sqrt{\frac{6}{d_{l}}}) , \sqrt{\frac{6}{d_{l}}})$ or $ \mathcal{N}(w;0,\sqrt{\frac{2}{d_{l}}})$

一般には$d_{l-1}$をfan-inに置き換えます。

効率的な勾配計算

ニューラルネットにおける入力データと推論結果を式で記述すると、以下のような巨大な合成関数になります。

$$\hat{y} = f(w^{(L)}f(w^{(L-1)}f(\dots w^{(l+1)}f(w^{(l)}f(\dots f(x)))\dots))$$

したがって、この関数を直接$w^{(l)}$で微分しようとすると途轍もない手間になってしまいます。

そこで、この面倒な計算をChain-Rule(合成関数の微分法)を用いて予め手で行っておき、出力層に近い部分から重みベクトルの勾配計算を行おうというのが「誤差逆伝播法(back propagation)」です。

今回はひとまず結果だけの紹介にとどめ、詳細の解説は別の機会に行おうと思います。

チェインルールより、

$$\frac{ \partial E}{ \partial w^{(l)}} = \frac{ \partial E}{ \partial u^{(l)}} \frac{ \partial u^{(l)}}{ \partial w^{(l)}} = \frac{\partial E}{\partial u^{(l)}} z^{(l-1)} \equiv \delta^{(l)} z^{(l-1)}$$

ここに$\delta^{(l)}_j = \frac{\partial E}{\partial u^{(l)}_j}$です。

また、$\delta^{(l)}_j$は次のようにして帰納的に求めることができます。

$$\delta^{(l)}_j =\sum_{k=0}^{d_{l+1}-1} \delta_k^{(l+1)}w_{kj}^{(l+1)}f^{‘}(u_j^{(l)})$$

以上をまとめると、以下の公式を得ます。

$$\frac{\partial E}{\partial w^{(l)}} = \delta^{(l)}(\textbf{z}^{(l-1)})^T$$

pythonによる実装

実装内容

主な内容は以下の三つです。
これに加えて、損失関数と正答率描写用の関数、学習用の関数、予測に用いる関数を実装しました。

  1. 誤差逆伝播法
  2. 勾配降下法
  3. gradient clipping

ソースコード

import numpy as np
from matplotlib import pyplot as plt
"""
2値分類用の、入力層ー中間層ー出力層の三層からなるパーセプトロンです。
"""

class MLP:

    def __init__(self, input_layer_size, hidden_layer_size, output_layer_size, initialize_storategy = 'random'):
        init_weight = 0.01 # 重みの初期値のばらつき
        
        # 今回はランダムに初期化
        self.W1 = init_weight*np.random.randn(input_layer_size, hidden_layer_size)
        self.b1 = np.zeros(hidden_layer_size)
        self.a1 = np.array([])
        self.z1 = np.array([])

        self.W2 = init_weight*np.random.randn(hidden_layer_size, output_layer_size)
        self.b2 = np.zeros(output_layer_size)
        self.a2 = np.array([])
        self.z2 = np.array([])

        self.loss_ = -1
        self.loss_curve_ = []
        self.accuracy_curve_ = []

 
    def fit(self, iters_num, learning_rate, X, Y, algorithm = 'gradientDescent', threshould = 0.5):
        for i in range(iters_num):
            Y_pred = self.forward_propagation(X) # 予測
            grad = self.backward_propagation(X, Y_pred, Y) # 勾配計算

            # 重みベクトル更新
            if algorithm == 'gradientDescent':
                self.gradient_descent(learning_rate, grad)
            elif algorithm == 'gradientClipping':
                self.gradient_clipping(learning_rate, grad, threshould)
            else:
                print(algorithm + 'is not supported.')
                return None

            self.accuracy_curve_.append(self.accuracy(X, Y))
            self.loss_curve_.append(self.loss(X, Y))
    
    def predict(self, X):
        return np.argmax(self.forward_propagation(X),axis=1)
    
    def predict_accuracy(self, X, y_real):
        y_real = np.argmax(y_real, axis = 1)
        y_pred = self.predict(X)
        return np.sum(y_pred == y_real)/float(X.shape[0])

# core algorithm
    def forward_propagation(self, x):
        self.a1 = np.dot(x, self.W1)+self.b1
        self.z1 = self.sigmoid(self.a1)
        self.a2 = np.dot(self.z1, self.W2)+self.b2
        self.z2 = self.sigmoid(self.a2)
        return self.z2

    def backward_propagation(self, x, y_pred, y_real):
        grads = Gradient()
        batch_num = x.shape[0]
        
        delta3 = (y_pred-y_real)/batch_num
        grads.W2 = np.dot(self.z1.T, delta3)
        grads.b2 = np.sum(delta3, axis = 0)
        
        delta2 = np.dot(delta3, self.W2.T)
        delta2 *= self.sigmoid_gradient(x = self.a1)
        grads.W1 = np.dot(x.T, delta2)
        grads.b1 = np.sum(delta2, axis = 0)
        
        return grads

# gradient descent
    def gradient_descent(self, learning_rate, grad):
        self.W1 -= learning_rate * grad.W1 
        self.b1 -= learning_rate * grad.b1 
        self.W2 -= learning_rate * grad.W2 
        self.b2 -= learning_rate * grad.b2 
    
    def gradient_clipping(self, learning_rate, grad, threshould):
        norm1 =  np.linalg.norm(grad.W1, ord=2)
        norm2 =  np.linalg.norm(grad.W2, ord=2)
        learning_rate1 = learning_rate
        learning_rate2 = learning_rate

        if norm1 > threshould:
            learning_rate1 *= threshould/norm1
        
        if norm2 > threshould:
            learning_rate2 *= threshould/norm2
        self.W1 -= learning_rate1 * grad.W1 
        self.b1 -= learning_rate1 * grad.b1 
        self.W2 -= learning_rate2 * grad.W2 
        self.b2 -= learning_rate2 * grad.b2        


# calcurate loss and accuracy
    def loss(self, x, y_real):
        self.loss_ = self.sigmoid_cross_entropy(self.forward_propagation(x), y_real)
        return self.loss_
    
    def accuracy(self, x, y_real):
        y_real = np.argmax(y_real, axis = 1)
        y_pred = np.argmax(self.forward_propagation(x), axis = 1)
        return np.sum(y_pred == y_real)/float(x.shape[0])


# plot loss and accuracy curve
    def plot_loss(self):
        plt.title("Loss Curve")
        plt.plot(self.loss_curve_)
        plt.xlabel("Iteration")
        plt.ylabel("Loss")
        plt.grid()
        plt.show()
    
    def plot_accuracy(self):
        plt.title("Accuracy Curve")
        plt.plot(self.accuracy_curve_)
        plt.xlabel("Iteration")
        plt.ylabel("Accuracy")
        plt.grid()
        plt.show()  

# define useful functions
    def sigmoid(self, x):
        return 1/(1+np.exp(-x))
    
    def sigmoid_gradient(self, x):
        return self.sigmoid(x)*(1-self.sigmoid(x))
    
    def sigmoid_cross_entropy(self, y_pred, y_real):
        if y_real.ndim == 1:
            y_real = np.reshape(1, y_real.size)
            y_pred = np.reshape(1, y_pred.size)
        
        if y_real.size == y_pred.size:
            y_real = y_real.argmax(axis = 1)
        
        batch_size = y_pred.shape[0]
        return -np.sum(np.log(y_pred[np.arange(batch_size), y_real] + 1e-8)) / batch_size
    


class Gradient:
    """
    勾配ベクトルを管理するためのクラス
    """
    def __init__(self):
        self.W1 = np.array([])
        self.b1 = np.array([])
        self.W2 = np.array([])
        self.b2 = np.array([])

この記事は役に立ちましたか?

もし参考になりましたら、下記のボタンで教えてください。

関連記事

コメント

この記事へのコメントはありません。

CAPTCHA