多層パーセプトロンによる関数近似 - 人工知能に関する断創録

パターン認識と機械学習（PRML）まとめ（2010/8/29）の続きです。以下、つづくかも？になってましたが、2014年はDeep Learningを勉強しよう（2014/1/4）と思っているので、関連するニューラルネットワーク関係の実験結果をもう少し追記します。

今回は、PRMLの5章ニューラルネットワークの中から図5.3にある多層パーセプトロンによる関数近似をPythonで実装してみました。

上の図にあるようにxを入れたときにsin(x)の近似値を出力するようなニューラルネットワークを学習します。もっと詳しく言うと (x, sin(x)) のたくさんの組（訓練データ）を教師データとして用いて、xを入れたときにsin(x)の近似値を出力するようにニューラルネットワークの重みを更新します。

多層パーセプトロン（Multilayer perceptron）

$(5.7') \hspace{50pt} \displaystyle y_k(\bf{x}, \bf{w}) = f \Biggl( \sum_{j=1}^M w_{kj}^{(2)} h \biggl( \sum_{i=1}^D w_{ji}^{(1)} x_i + w_{j0}^{(1)} \biggr) + w_{k0}^{(2)} \Biggr)$

p.228にあるようにバイアスパラメータは、入力変数にx0=1を追加することで重みパラメータの中に含められるので

$(5.9') \hspace{50pt} \displaystyle y_k(\bf{x}, \bf{w}) = f \Biggl( \sum_{j=0}^M w_{kj}^{(2)} h \biggl( \sum_{i=0}^D w_{ji}^{(1)} x_i \biggr) \Biggr)$

のようにすっきりした式になります。出力ユニット活性化関数はf()、隠れユニット活性化関数はh()で表しています。PRMLの式(5.7)と式(5.9)は出力ユニット活性化関数にシグモイド関数を用いていますがここではより一般的なf()で書いています。今回の例では、図5.3に書いてあるように出力ユニット活性化関数には線形関数、隠れユニット活性化関数にはtanh()を用います。

Dは入力層のユニット数（バイアス除く）、Mは隠れ層のユニット数（バイアス除く）です。今回は、任意の数字xを入力するのでD=1です。隠れ層のユニット数は適当にM=3としています。数字f(x)が1個出てくればいいので出力層のユニット数は1です。よってバイアスユニット（灰色の丸）を含めて下の図のような多層パーセプトロンを構成しました。 $x_0$ は必ず1を入れます。

最初、バイアスユニットを入れ忘れていたため正しい結果が得られず、ひどく悩みました。バイアスユニットは式 (5.7) の定数項を表現しているためないとダメなんですね・・・定数項も重要なのだ。

後述するプログラムは中間層のバイアス項の出力z0はそのままx0とx1の重み付け和で計算してしまっていました。上の式からするに順伝播における中間層のバイアスユニットの出力も常に1.0にしておくのが正しいのかもしれません。なのでzを計算し終わった後に z[0] = 1.0 という式が必要かも。まあそうしなくても収束してしまいますが・・・中間層を多めにしているからかな？
（2014/1/27 追記）

誤差逆伝播法（Backpropagation）

重みの学習法には、全学習データの二乗誤差の和

$(5.44) \hspace{50pt} \displaystyle E(\bf{w}) = \sum_{n=1}^N E_n(\bf{w})$
$(5.61) \hspace{50pt} \displaystyle E_n = \frac{1}{2} \sum_{k=1}^K (y_k - t_k)^2$

を最小化する誤差逆伝播法を用います。誤差逆伝播法のアルゴリズムは、

(1) 順伝播

入力ベクトル $x_n$ をネットワークに入れ、すべての隠れユニットの出力 (5.63) と出力ユニットの出力 (5.64) を求める。

$(5.63) \hspace{50pt} \displaystyle z_j = \tanh \biggl( \sum_{i=0}^D w_{ji}^{(1)} x_i \biggr)$
$(5.64) \hspace{50pt} \displaystyle y_k = \sum_{j=0}^M w_{kj}^{(2)} z_j$

(2) 誤差の計算

(5.65) を用いてすべての出力ユニットの誤差 $\delta_k$ を求める。

$(5.65) \hspace{50pt} \delta_k = y_k - t_k$

(3) 逆伝播

(5.66) を用いて誤差を逆伝播させ、ネットワークのすべての隠れユニットの誤差 $\delta_j$ を求める。 $(1 - z_j^2)$ は、tanh()の微分で隠れ層の活性化関数にtanh()を使っているため。sigmoid()など他の活性化関数を使う場合は別の式になるので注意。

$(5.66) \hspace{50pt} \displaystyle \delta_j = (1 - z_j^2) \sum_{k=1}^K w_{kj} \delta_k$

(4) 誤差の微分

(5.67) を用いて第1層と第2層の重みに関する微分を求める。

$(5.67) \hspace{50pt} \displaystyle \frac{\partial E_n}{\partial w_{ji}^{(1)}} = \delta_j x_i, \hspace{30pt} \frac{\partial E_n}{\partial w_{kj}^{(2)}} = \delta_k z_j$

(5) 重みの更新

(5.43) を用いて第1層と第2層の重みをそれぞれ更新する。 $\nabla E_n(\bf{w})$ の部分は上の式 (5.67) です。また、 $\eta$ は学習率。

$(5.43) \hspace{50pt} \bf{w}^{(\tau + 1)} = \bf{w}^{(\tau)} - \eta \nabla E_n (\bf{w}^{(\tau)})$

以上の過程をすべてのデータに対して繰り返し行い、式 (5.44) の全学習データの二乗誤差の和が十分小さくなったら学習を終了する。

Pythonによる実装

というわけでさっそくプログラムを書いて本当に学習できるか確かめてみます。少し効率は悪いですが、上の式に忠実に実装してみます。

#coding:utf-8
import numpy as np
import pylab

# 訓練データ数
N = 50

# 学習率
ETA = 0.1

# ループ回数
NUM_LOOP = 1000

# 入力層のユニット数（バイアス含む）
NUM_INPUT = 2
# 隠れ層のユニット数（バイアス含む）
NUM_HIDDEN = 4
# 出力層のユニット数
NUM_OUTPUT = 1

def heviside(x):
    """Heviside関数"""
    return 0.5 * (np.sign(x) + 1)

def sum_of_squares_error(xlist, tlist, w1, w2):
    """二乗誤差和を計算する"""
    error = 0.0
    for n in range(N):
        z = np.zeros(NUM_HIDDEN)
        y = np.zeros(NUM_OUTPUT)
        # バイアスの1を先頭に挿入
        x = np.insert(xlist[n], 0, 1)
        # 順伝播で出力を計算
        for j in range(NUM_HIDDEN):
            a = np.zeros(NUM_HIDDEN)
            for i in range(NUM_INPUT):
                a[j] += w1[j, i] * x[i]
            z[j] = np.tanh(a[j])

        for k in range(NUM_OUTPUT):
            for j in range(NUM_HIDDEN):
                y[k] += w2[k, j] * z[j]
        # 二乗誤差を計算
        for k in range(NUM_OUTPUT):
            error += 0.5 * (y[k] - tlist[n, k]) * (y[k] - tlist[n, k])
    return error

def output(x, w1, w2):
    """xを入力したときのニューラルネットワークの出力を計算
    隠れユニットの出力も一緒に返す"""
    # 配列に変換して先頭にバイアスの1を挿入
    x = np.insert(x, 0, 1)
    z = np.zeros(NUM_HIDDEN)
    y = np.zeros(NUM_OUTPUT)
    # 順伝播で出力を計算
    for j in range(NUM_HIDDEN):
        a = np.zeros(NUM_HIDDEN)
        for i in range(NUM_INPUT):
            a[j] += w1[j, i] * x[i]
        z[j] = np.tanh(a[j])
    for k in range(NUM_OUTPUT):
        for j in range(NUM_HIDDEN):
            y[k] += w2[k, j] * z[j]
    return y, z

if __name__ == "__main__":
    # 訓練データ作成
    xlist = np.linspace(-1, 1, N).reshape(N, 1)
    tlist = xlist * xlist    # x^2
#    tlist = np.sin(xlist)    # sin(x)
#    tlist = np.abs(xlist)    # |x|
#    tlist = heviside(xlist)  # H(x)

    # 重みをランダムに初期化
    w1 = np.random.random((NUM_HIDDEN, NUM_INPUT))
    w2 = np.random.random((NUM_OUTPUT, NUM_HIDDEN))

    print "学習前の二乗誤差:", sum_of_squares_error(xlist, tlist, w1, w2)

    # 収束するまで十分なループを回す
    # 二乗誤差和がepsilon以下になったら終了でもOK
    for loop in range(NUM_LOOP):
        # 訓練データすべてを使って重みを訓練する
        for n in range(len(xlist)):
            # 隠れ層と出力層の出力配列を確保
            z = np.zeros(NUM_HIDDEN)
            y = np.zeros(NUM_OUTPUT)

            # 誤差（delta）の配列を確保
            d1 = np.zeros(NUM_HIDDEN)
            d2 = np.zeros(NUM_OUTPUT)

            # 訓練データにバイアスの1を先頭に挿入
            x = np.array([xlist[n]])
            x = np.insert(x, 0, 1)

            # (1) 順伝播により隠れ層の出力を計算
            for j in range(NUM_HIDDEN):
                # 入力層にはバイアスの1が先頭に入るので注意
                a = np.zeros(NUM_HIDDEN)
                for i in range(NUM_INPUT):
                    a[j] += w1[j, i] * x[i]
                z[j] = np.tanh(a[j])

            # (1) 順伝播により出力層の出力を計算
            for k in range(NUM_OUTPUT):
                for j in range(NUM_HIDDEN):
                    y[k] += w2[k,j] * z[j]

            # (2) 出力層の誤差を評価
            for k in range(NUM_OUTPUT):
                d2[k] = y[k] - tlist[n, k]

            # (3) 出力層の誤差を逆伝播させ、隠れ層の誤差を計算
            for j in range(NUM_HIDDEN):
                temp = 0.0
                for k in range(NUM_OUTPUT):
                    temp += w2[k, j] * d2[k]
                d1[j] = (1 - z[j] * z[j]) * temp

            # (4) + (5) 第1層の重みを更新
            for j in range(NUM_HIDDEN):
                for i in range(NUM_INPUT):
                    w1[j, i] -= ETA * d1[j] * x[i]
            
            # (4) + (5) 第2層の重みを更新
            for k in range(NUM_OUTPUT):
                for j in range(NUM_HIDDEN):
                    w2[k, j] -= ETA * d2[k] * z[j]

        # 全データで重み更新後の二乗誤差和を出力
        print loop, sum_of_squares_error(xlist, tlist, w1, w2)

    # 学習後の重みを出力
    print "w1:", w1
    print "w2:", w2

    # 訓練データの入力に対する隠れユニットと出力ユニットの出力を計算
    ylist = np.zeros((N, NUM_OUTPUT))
    zlist = np.zeros((N, NUM_HIDDEN))
    for n in range(N):
        ylist[n], zlist[n] = output(xlist[n], w1, w2)

    # 訓練データを青丸で表示
    # plot()には縦ベクトルを渡してもOK
    pylab.plot(xlist, tlist, 'bo')

    # 訓練データの各xに対するニューラルネットの出力を赤線で表示
    pylab.plot(xlist, ylist, 'r-')

    # 隠れユニットの出力を点線で表示
    for i in range(NUM_HIDDEN):
        pylab.plot(xlist, zlist[:,i], 'k--')

    pylab.ylim([0, 1])
    pylab.show()

ループを回すところなどはnumpy.dot()で2つのベクトルの内積計算に置き換えた方がより効率的です。例えば、順伝播のところは

            # (1) 順伝播により隠れ層の出力を計算
            a = np.dot(w1, x)
            z = np.tanh(a)

            # (1) 順伝播により出力層の出力を計算
            y = np.dot(w2, z)

            # (2) 出力層の誤差を評価
            d2 = y - tlist[n,:]

と行列演算しても同じ結果になります。

順伝播、逆伝播含めほとんどのループは行列演算に置き換えることができます。
http://rolisz.ro/2013/04/18/neural-networks-in-python/
のコードはPRMLとweightsの格納方法や添字が違いますが参考になります。
多層パーセプトロンで手書き数字認識（2014/2/1）
http://aidiary.hatenablog.com/entry/20140201/1391218771
では行列演算を駆使して再実装しました。

PRMLの図5.3と同様の関数で試してみた結果です。関数によってループ回数や表示範囲を調整しています。青い点が訓練データで赤線がニューラルネットの出力です。黒い点線は隠れユニットの出力です。黒い点線を第2層の重み（w2）で重み付けして和を取った値が赤い線ですね。

f:id:aidiary:20140122212500p:plain
f:id:aidiary:20140122212505p:plain
f:id:aidiary:20140122212509p:plain
f:id:aidiary:20140122212513p:plain

一番最後のヘヴィサイドステップ関数は収束が遅かったです。結果もいまいちですね・・・また、重みの初期値はランダムにしていますが、初期値によっては収束が非常に遅い場合もありました。何回か試してみるとよいかも。あと隠れ層のユニット数や学習率を変えることで、学習結果がどのように変わるか分析するのも面白いかも。

次回は、上記の多層パーセプトロンの収束する様子がよくわかるようにアニメーションで表示（2014/1/23）してみたいと思います。