Deep Learning のための Multi Layer Perceptron (数学的基礎から学ぶ Deep Learning with Python; MPS Yokohama Deep Learning Series)

目的

MLP (Multi Layer Perceptron) を Python3 で Numpy と Scipy のみを使って作成する。また、実際の例として手書き数字データベース MNIST を用いて、手書き数字画像のクラス分類を行う MLP の構築を行う。

環境

ソフトウェア	バージョン
Python	3.5
Numpy	1.10
Scipy	0.16
Matplotlib	1.5
tqdm	4.7
python-mnist	0.3

注意

本記事は、2016/7/23 に行われる MPS (Morning Project Samurai) Yokohama のイベント 数学的基礎から学ぶ Deep Learning (with Python) Vol. 10 のために書かれた記事である。質問や指摘等あれば頂けると幸いである。

MLP (Multi Layer Perceptron)

まず、下図に示すような単純な MLP を考える。実際には layer $i + 1$ と layer $i$ の node は全て結合されている (fully connected) が、見易さのため結線の一部しか描いていない。

ここで、下記のように記号の定義をしておく。

• $s_{ij}$: layer $i$ の $j$ 番目の node の状態
• $w_{ikj}$: layer $i$ の $k$ 番目の node に入力される layer $i-1$ の $j$ 番目の出力にかかる重み
• $f_{i}$: layer $i$ の AF (Activation Function)
• $y_{ij}$: layer $i$ の $j$ 番目の node の出力
• $E$: EF (Error Function)

FP (Forward Propagation)

layer $i-1$ から layer $i$ への FP は次のようになる。

y_{i} = f_{i}(W_{i}y_{i-1} + b_{i-1})

ここで $y_{i} = (y_{ij})$, $W_{i}=(w_{ikj})$, $b_{i}=(b_{ij})$ とする。

3入力1出力を持つ2層の MLP の FP の Python コードは次のようになる。

import numpy as np


# activation function
def sigmoid(s):
    return 1/(1 + np.exp(-s))


# x: Input vector, W: Weight matrix, b: Bias, y: Output vector
x = np.random.randn(3, 1)
W = np.random.randn(1, 3)
b = np.random.randn(1, 1)
y = sigmoid(W @ x + b)

上記オペレーションをメソッドに持つクラス Layer を定義すれば、MLP を簡単に作成できるようになる。下記は入力層と出力層に加え、1つの中間層を持ったネットワークの例である。

from scipy import misc
import numpy as np


class Layer:
    def __init__(self, W, b, f):
        self._W = W
        self._b = b
        self._f = f

    def propagate_forward(self, x):
        return self._f(self._W @ x + self._b)


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


if __name__ == '__main__':
    # Input vector (Layer 0)
    x = np.random.randn(3, 1)

    # Middle layer (Layer 1)
    W1 = np.random.randn(3, 3)
    b1 = np.random.randn(3, 1)
    layer1 = Layer(W1, b1, sigmoid)

    # Output layer (Layer 2)
    W2 = np.random.randn(3, 3)
    b2 = np.random.randn(3 ,1)
    layer2 = Layer(W2, b2, sigmoid)

    # FP
    y1 = layer1.propagate_forward(x)
    y2 = layer2.propagate_forward(y1)

    print(y2)

ここまでで、MLP の構成と、その MLP のある入力に対する出力を得られるようになった。以下は、下図のような手書き数字の9を入力として与えた例である。

from scipy import misc
import numpy as np
from matplotlib import pyplot as plt


class Layer:
    def __init__(self, W, b, f):
        self._W = W
        self._b = b
        self._f = f

    def propagate_forward(self, x):
        return self._f(self._W @ x + self._b)


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


if __name__ == '__main__':
    # Input vector (Layer 0)
    x = misc.imread('img/9.png', flatten=True).flatten()
    x = x.reshape((len(x), 1))

    # Middle layer (Layer 1)
    n_output_1 = len(x)
    W1 = np.random.randn(n_output_1, len(x))
    b1 = np.random.randn(n_output_1, 1)
    layer1 = Layer(W1, b1, sigmoid)

    # Output layer (Layer 2)
    n_output_2 = 10
    W2 = np.random.randn(n_output_2, n_output_1)
    b2 = np.random.randn(n_output_2, 1)
    layer2 = Layer(W2, b2, sigmoid)

    # FP
    y1 = layer1.propagate_forward(x)
    y2 = layer2.propagate_forward(y1)

    hist_W1, bins_W1 = np.histogram(W1.flatten())
    hist_W2, bins_W2 = np.histogram(W2.flatten())

    # Draw bar chart
    index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    plt.title('Prediction')
    plt.bar(index , y2.flatten(), align="center")
    plt.xticks(index, index)
    plt.show()

Activation function に sigmoid 関数を使用しているので、出力層にある 10 個の node の出力は 0 から 1 の間をとる。出力層の $i$ 番目の node の出力を入力画像が数字 $i$ である確率を表しているとすると、上記結果の場合、2, 6 または 7 である確率が最も高いと予測しており、役に立たない。

EF (Error Function)

EF は MLP の状態が理想的な状態からどの程度かけ離れているかを測るための指標である。例えば、上述の手書き画像認識の場合、正解は 9 であり出力として下記を得るのが理想である。

y_{2} = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0)^{T}

よく用いられる誤差関数の1つとして SE (Squared Error) があり、下記のように表される。

E = \frac{1}{2}\sum_{j} (t_{j} - y_{ij})^{2} = \frac{1}{2}(t-y_{i})^{T}(t-y_{i})

手書きの 9 を入力としたときの SE を求めるコードを下記に示す。

from scipy import misc
import numpy as np


class Layer:
    def __init__(self, W, b, f):
        self._W = W
        self._b = b
        self._f = f

    def propagate_forward(self, x):
        return self._f(self._W @ x + self._b)


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


def se(t, y):
    return ((t - y).T @ (t-y)).flatten()[0]/2


if __name__ == '__main__':
    # Input vector (Layer 0)
    x = misc.imread('img/9.png', flatten=True).flatten()
    x = x.reshape((len(x), 1))

    # Middle layer (Layer 1)
    n_output_1 = len(x)
    W1 = np.random.randn(n_output_1, len(x))
    b1 = np.random.randn(n_output_1, 1)
    layer1 = Layer(W1, b1, sigmoid)

    # Output layer (Layer 2)
    n_output_2 = 10
    W2 = np.random.randn(n_output_2, n_output_1)
    b2 = np.random.randn(n_output_2, 1)
    layer2 = Layer(W2, b2, sigmoid)

    # FP
    y1 = layer1.propagate_forward(x)
    y2 = layer2.propagate_forward(y1)

    # Calculate SE
    t = np.zeros(shape=(10, 1))
    t[9, 0] = 1.0
    e = se(t, y2)
    print(e)

BP (Backward Propagation) と学習

BP は、情報をネットワークの逆向きに伝搬させてゆくことである。BP を用いた学習では、学習に必要な情報を BP を用いて逆向きに伝搬し、その情報を用いて実際に学習を行うという2ステップが必要となる。

MLP の学習は、EF の gradient を学習させたいパラメータを domain とする空間で求め、その反対方向へパラメータを移動させることが現在のパラメータの近傍において最も EF の値を減少させる、という事実を用いて行う。BP は、この gradient の計算を効率的に行うために必要な情報の伝搬を行うために用いられる。

出力層の学習

まず最初に、3層の MLP の出力層の学習を考える。

調整したいパラメータは重み $w_{2ji}$ であるので、それで EF $E$ の偏微分をすると次式が得られる。

\frac{\partial{E}}{\partial{w_{2ji}}} = 
\frac{\partial{E}}{\partial{y_{2j}}}
\frac{\partial{y_{2j}}}{\partial{s_{2j}}}
\frac{\partial{s_{2j}}}{\partial{w_{2ji}}}

これは、$w_{2ji}$ を変化させたときに $E$ がどのように変化するかを調べるには、$w_{2ji}$ を変化させた時に $s_{2j}$ がどのように変化するかを調べ、次に $s_{2j}$ を変化させた時に $y_{2j}$ がどのように変化するかを調べ、最後に $y_{2j}$ を変化させた時に $E$ がどのように変化するかを調べれば良い、という関係を表しており、MLP の図と照らし合わせて見れば直感的なものとなっている。

ここで、次のようにおく。

\frac{\partial{E}}{\partial{y_{2j}}}\frac{\partial{y_{2j}}}{\partial{s_{2j}}} = \delta_{2j}

また、$s_{2j} = \sum_{i} w_{2ji}y_{1i}$ より次式が得られる。

\frac{\partial{s_{2j}}}{\partial{w_{2ji}}} = y_{1i}

これまでの式を合わせて次式を得る。

\frac{\partial{E}}{\partial{w_{2ji}}} = \delta_{2j}y_{1i}

学習は次の更新式で行うことができる。ここで式中の矢印は、左辺の変数を右辺で置き換えることを意味する。また、 $\epsilon$ は学習率と呼ばれるハイパーパラメータである。この値を調整することにより、学習の速度を速めたり、遅くしたりすることができる。

w_{2ji} \leftarrow w_{2ji} - \epsilon \delta_{2j}y_{1i}

ここで、次式のようにおき、

\delta_{2} = (\delta_{2j})

これまでの議論を行列の形で整理すると次式のようになる。

\Delta W_{2} = (\frac{\partial{E}}{\partial{w_{2ji}}}) 
= \delta_{2}y_{1}^{T}

W_{2} \leftarrow W_{2} - \epsilon \Delta W_{2}

Bias 項についても同様に考えれば、次式が得られる。

\frac{\partial{E}}{\partial{b_{2j}}} = \delta_{2j}

よって bias 項の更新式は次式のようになる。

b_{2} \leftarrow b_{2} - \epsilon \delta_{2}

具体的に AF として sigmoid 関数を、EF として SE を用いた例を示す。

\frac{\partial{E}}{\partial{y_{2j}}} = -(t_{j}-y_{2j})

\frac{\partial{y_{2j}}}{\partial{s_{2j}}} = y_{2j}(1-y_{2j})

よって、次式が得られる。

\delta_{2j} = -(t_{j}-y_{2j})y_{2j}(1-y_{2j})

この $\delta_{2} = (\delta_{2j})$ を用いれば、重みと Bias 両方のパラメータの学習を行うことができる。

下記は、3層の MLP に手書き文字 9 の認識を出力層のみで学習させるコードである。

from scipy import misc
import numpy as np
from matplotlib import pyplot as plt


class Layer:
    def __init__(self, W, b, f):
        self._W = W
        self._b = b
        self._f = f

    def propagate_forward(self, x):
        return self._f(self._W @ x + self._b)


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


def d_sigmoid(y):
    return y * (1 - y)


def se(t, y):
    return ((t - y).T @ (t-y)).flatten()[0]/2


def d_se(t, y):
    return -(t - y)


if __name__ == '__main__':
    # Input vector (Layer 0)
    x = misc.imread('img/9.png', flatten=True).flatten()
    x = x.reshape((len(x), 1))

    # Middle layer (Layer 1)
    n_output_1 = len(x)
    W1 = np.random.randn(n_output_1, len(x))
    b1 = np.random.randn(n_output_1, 1)
    layer1 = Layer(W1, b1, sigmoid)

    # Output layer (Layer 2)
    n_output_2 = 10
    W2 = np.random.randn(n_output_2, n_output_1)
    b2 = np.random.randn(n_output_2, 1)
    layer2 = Layer(W2, b2, sigmoid)

    # Training datum
    t = np.zeros(shape=(10, 1))
    t[9, 0] = 1.0

    # FP, BP and learning
    epsilon = 0.1
    se_history = []
    delta2_history = []
    for i in range(0, 100):
        # FP
        y1 = layer1.propagate_forward(x)
        y2 = layer2.propagate_forward(y1)

        # Calculate and store SE
        se_history.append(se(t, y2))

        # BP and learning
        delta2 = d_se(t, y2) * d_sigmoid(y2)
        Delta_W2 = delta2 @ y1.T
        layer2._W -= epsilon * Delta_W2
        layer2._b -= epsilon * delta2

        # Store delta2 history
        delta2_history.append(np.linalg.norm(delta2))

    y1 = layer1.propagate_forward(x)
    y2 = layer2.propagate_forward(y1)

    # Draw SE history
    plt.figure()
    plt.title('SE History')
    plt.plot(range(len(se_history)), se_history)

    # Draw delta2 history
    plt.figure()
    plt.title('delta2 history')
    plt.plot(range(len(delta2_history)), delta2_history)

    # Draw bar chart
    index = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
    plt.figure()
    plt.title('Prediction')
    plt.bar(index, y2.flatten(), align="center")
    plt.xticks(index, index)

    plt.show()

下図は、100回学習させた時の SE の推移を表したものである。学習の初期段階で大きく SE が下がっていることが見て取れる。また、5回目くらい以降はほとんど学習の効果が出ていないこともわかる。

その原因を確かめるために、$\delta_{2}$ の大きさのグラフを下記に示す。5回目くらい以降から、その大きさがほぼ $0$ になっていることがわかる。

この時の MLP 出力は次のようになっている。

中間層の学習

次に、3層 MLP の中間層の学習を考える。下にもう一度 3 層 MLP の図を示す。

出力層の学習と同様に、まず学習させたい重み $w_{1ih}$ で EF を偏微分する。ここで、例えば $w_{100}$ を変化させた場合、この変化は $s_{10}$ に影響を及ぼし、次いで $y_{10}$ に影響を及ぼし、これは出力層の全ての node $s_{2j}$ に影響を及ぼし、それらの出力 $y_{2j}$ に影響を及ぼし、最終的にそれらを通じて EF の値に影響を及ぼす。これを考慮して偏微分の式を展開し、整理してゆく。

\begin{eqnarray}
\frac{\partial{E}}{\partial{w_{1ih}}} 
&=& \sum_{j}\frac{\partial{E}}{\partial{y_{2j}}}
\frac{\partial{y_{2j}}}{\partial{s_{2j}}}
\frac{\partial{s_{2j}}}{\partial{y_{1i}}}
\frac{\partial{y_{1i}}}{\partial{s_{1i}}}
\frac{\partial{s_{1i}}}{\partial{w_{1ih}}}\\

&=& (\sum_{j}\frac{\partial{E}}{\partial{y_{2j}}}
\frac{\partial{y_{2j}}}{\partial{s_{2j}}}
\frac{\partial{s_{2j}}}{\partial{y_{1i}}})
\frac{\partial{y_{1i}}}{\partial{s_{1i}}}
\frac{\partial{s_{1i}}}{\partial{w_{1ih}}}\\

&=&(\sum_{j}\delta_{2j} w_{2ji})
\frac{\partial{y_{1i}}}{\partial{s_{1i}}}y_{0h}\\

&=&\delta_{1i}y_{0h}
\end{eqnarray}

ここで、$\delta_{1i}$ を次のようにおいた。

\delta_{1i} = (\sum_{j}\delta_{2j} w_{2ji})
\frac{\partial{y_{1i}}}{\partial{s_{1i}}}

上の式変形の最後で得た式は、出力層の重みの学習を考えた時に導いた式と形が全く同じであることに注意する。それらの式を再掲する。

\begin{eqnarray}
\frac{\partial{E}}{\partial{w_{2ji}}} &=& \delta_{2j}y_{1i}\\
\frac{\partial{E}}{\partial{w_{1ih}}} &=& \delta_{1i}y_{0h}
\end{eqnarray}

上が出力層の学習で用いたもの、下が今回の中間層の学習で用いるものである。

$\delta_{1i}$ は出力層の学習のために求めた $\delta_{2j}$ をもとに作られていることにも注意する。つまり、出力層で求めた $\delta_{2j}$ さえ分かっていれば (出力層から伝搬されてきていれば)、複雑な計算をすることなく中間層の $\delta_{1i}$ を求めることができ、それをもとに中間層の学習ができることになる。MLP の学習では、この $\delta$ を誤差と呼び、これを伝搬することを BP と呼んでいる。

中間層の重みの学習を行列表現を用いて整理しておく。

\delta_{1} = W_{2}^{T} \circ \frac{\partial{y_{1}}}{\partial{s_{1}}}

ここで、$\circ$ はアダマール積を表す。アダマール積は、行列の要素ごとに掛け算を行うものである。

\Delta W_{1} = \delta_{1} y_{0}^T

ここで $y_{0}$ は入力層の各 node の出力で構成されるベクトルである。学習による重みの更新式は次のようになる。

W_{1} \leftarrow W_{1} - \epsilon \Delta W_{1}

Bias 項については、是非各自導出を試みて頂きたい。

手書き数字画像のクラス分類を行う MLP の作成

ここまでの議論をもとに MNIST と呼ばれる手書き数字のデータベースを用いて 3 層 MLP に手書き数字のクラス分類を学習させるコードを次に示す。これまでと同様、出力層には 10 個の node があり、$i$ 番目の node は与えられた画像に書かれた数字が $i$ である確率を表すとする。今回は最も高い確率を示した node の番号を最終的な認識結果として採用する。

MNIST は、60000 万個の手書き数字の画像をラベルとともに教師データとして提供しており、今回はこのデータの最初の 1000 個を学習に用い、その他をテストデータとして用いている。学習回数は上記データを 400 回繰り返し用い、合計 40 万回としてる。

import numpy as np
from matplotlib import pyplot as plt
from tqdm import tqdm


class Layer:
    def __init__(self, W, b, f):
        self._W = W
        self._b = b
        self._f = f

    def propagate_forward(self, x):
        return self._f(self._W @ x + self._b)


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


def d_sigmoid(y):
    return y * (1 - y)


def se(t, y):
    return ((t - y).T @ (t - y)).flatten()[0] / 2.0


def d_se(t, y):
    return -(t - y)


def ma(history, n):
    return np.array([0, ] * (n - 1) + [np.average(history[i - n: i]) for i in range(n, len(history) + 1)])


if __name__ == '__main__':
    from mnist import MNIST

    # Load MNIST dataset
    mndata = MNIST('./mnist')
    train_img, train_label = mndata.load_training()
    train_img = np.array(train_img, dtype=float)/255.0
    train_label = np.array(train_label, dtype=float)

    # Input vector (Layer 0)
    n_output_0 = len(train_img[0])

    # Middle layer (Layer 1)
    n_output_1 = 200
    W1 = np.random.randn(n_output_1, n_output_0)
    b1 = np.random.randn(n_output_1, 1)
    layer1 = Layer(W1, b1, sigmoid)

    # Output layer (Layer 2)
    n_output_2 = 10
    W2 = np.random.randn(n_output_2, n_output_1)
    b2 = np.random.randn(n_output_2, 1)
    layer2 = Layer(W2, b2, sigmoid)

    # FP, BP and learning
    epsilon = 0.15
    n_training_data = 1000
    se_history = []
    y1_history = []
    y2_history = []
    W1_history = []
    W2_history = []
    cpr_history = []
    for loop in range(400):
        for i in tqdm(range(n_training_data)):
            # Store W1 and W2 history
            W1_history.append(np.linalg.norm(layer1._W))
            W2_history.append(np.linalg.norm(layer2._W))

            # FP
            x = train_img[i].reshape(len(train_img[i]), 1)
            y1 = layer1.propagate_forward(x)
            y2 = layer2.propagate_forward(y1)

            # Store y1 and y2
            y1_history.append(y1)
            y2_history.append(y2)

            # Training datum
            t = np.zeros(shape=(10, 1))
            t[train_label[i], 0] = 1.0

            # Calculate and store SE
            se_history.append(se(t, y2))

            # BP
            delta2 = d_se(t, y2) * d_sigmoid(y2)
            delta1 = layer2._W.T @ delta2 * d_sigmoid(y1)

            # Learning
            Delta_W2 = delta2 @ y1.T
            layer2._W -= epsilon * Delta_W2
            layer2._b -= epsilon * delta2

            Delta_W1 = delta1 @ x.T
            layer1._W -= epsilon * Delta_W1
            layer1._b -= epsilon * delta1

        # FP to evaluate correct prediction rate
        n_correct_prediction = 0
        n_prediction = 0
        for _i in np.random.choice(np.arange(n_training_data, train_img.shape[0]), 100):
            _x = train_img[_i].reshape(len(train_img[_i]), 1)
            _y1 = layer1.propagate_forward(_x)
            _y2 = layer2.propagate_forward(_y1)

            n_prediction += 1
            if train_label[_i] == np.argmax(_y2):
                n_correct_prediction += 1
        cpr_history.append(n_correct_prediction/n_prediction)

    # Draw W1
    plt.figure()
    plt.title('W1 history')
    plt.plot(range(len(W1_history)), W1_history)

    # Draw W2
    plt.figure()
    plt.title('W2 history')
    plt.plot(range(len(W2_history)), W2_history)

    # Draw SE history and its moving average
    plt.figure()
    plt.title('SE History')
    plt.plot(range(len(se_history)), se_history, color='green')
    plt.plot(range(len(se_history)), ma(se_history, 100), color='red')

    # Draw CPR history
    plt.figure()
    plt.title('CPR')
    plt.plot(range(len(cpr_history)), cpr_history)

    plt.show()

下図は、MLP の SE の変化を学習回数を横軸に取って表したものである。緑は、1 つのデータを学習し、その後、他のデータが 1 つ与えられた時の SE の値を表している。赤は緑の系列の移動平均になっている。平均に使用しているデータは 100 個なので、100 回の学習で得られる SE の平均値の系列になっている。この図をみれば、学習を重ねるごとに、学習 1 回 1 回で得られる SE は大きく振動しているが、平均的に見ると小さくなっていることが見てとれる。

手書き数字の認識率の変化を学習回数を横軸に取って表した図を下に示す。ただし、予測精度測定は 1000 回の学習ごとに 1 度行っているので、実際の学習回数は横軸の値に 1000 を掛けたものになる。ここでも、予測精度は、移動平均をみれば学習を重ねるごとに向上していることが見て取れる。このことから、予測は 1 つの MLP だけで行うのではなく、複数の MLP を組み合わせて行うほうが良いことが予想される。

最後に

本記事で実装したのは基本的な MLP とその学習である。MPS Yokohama の 数学的基礎から学ぶ Deep Learning (with Python) では、次回から Deep Learning を行うようこれを改造してゆく。

本企画においては、有志の方々が素晴らしい資料を作ってくださっているので、そちらも是非参考にしていただきたい (順不同)。

関連資料	著者
MPS Yokohama 第8回までのまとめ	さいとうさん
Back Propagation の整理	n-ken さん
Python による微分のコード	takaneh さん
MPS Yoohama 第9回補講ノート	てらさきさん
合成関数の連鎖率	てさらきさん
1変数関数の連鎖率	てらさきさん
まとめと疑問点(作成中)	わたるさん

また SlideShare と YouTube 上にこれまでのイベントで使用したスライドとイベントの動画もアップしてあるので、ご興味ある方はご覧いただけると幸いである。

| 回 | スライド | 勉強会の動画 |
|:-------:|:----------:|:-------------:|:-------------:|
| 1 | Vol. 1 |Youtube|
| 2 | Vol. 2 | Youtube |
| 3 | Vol. 3 | Youtube |
| 4 | Vol. 4 |Youtube|
| 5 | Vol. 5 | YouTube |
| 6 | Vol. 6 | YouTube|
| 7 | Vol. 7 | Youtube |
| 8 | Vol. 8 | YouTube |
| 9 | Vol. 9 | Youtube |

謝辞

本企画は多くの方々の多大な協力によって成り立っています。毎回、会場やプロジェクター、WiFiといった設備を貸してくださり、ビデオ撮影とそのアップロードを行ってくださっています情報科学専門学校様、武藤先生、クラウド環境をご提供くださっていますサクラインターネット様、毎回のまとめを作ってくださり、また勉強会の運営に関して的確なアドバイスをくださる、さいとうさん、MPS Yokohama のリーダーとして本勉強会を引っ張ってくださっている上野さん、毎朝運営のお手伝いをしてくださる、わたるさん、資料を作成くださったり毎回質問や指摘をしてくださる、てらさきさん、ひらはらさん、n-ken さん、本当にありがとうございます。また、出会ってから2年間、ずっと MPS の活動を支えてくださっているまさきさん、本当にありがとうございます。この場を借りてお礼申し上げます。これからも、宜しくお願いいたします。

2016/07/22 Junya Kaneko