前回、特に数字の特徴を定義したわけでもないのに、97%まで認識率を向上させることができました。不思議でたまりません。

パラメータ

前回の実行では、パラメータは、適当に選択したため、今回は、パラメータを変更し、どのような違いがあるか確認していきます。データは、今まで通り、MNISTを利用します。
事前に少し実行し、今回試すパラメータを見直しました。実行時のパラメータを以下にまとめます。

パラメータ	前回	今回
階層数	3	3
1階層目のノード数	50	100
2階層目のノード数	100	50
重みの初期値	he_normal	he_normal
バイアスの初期値	0	0
正規化関数	min_max	min_max
活性化関数(中間層)	relu	relu
活性化監修(出力層)	softmax	softmax
損失関数	cross_entropy_error	cross_entropy_error
学習率	0.01	0.1
バッチサイズ	100	100
エポック数	100	50

パラメータ変更の組み合わせは無限にあるため、一項目ごと変更し試します。

実装変更

現在は、3階層固定で実装していました。今回、階層数を変更可能とするため実装を見直します。合わせて、所要時間と誤差を表示するようにしました。また、後で正解率、誤差のグラフが描けるように、エポックごとの正解率、誤差を保存します。

まずは、階層を自由に変更できるよう対応します。
各階層のノード数を配列で定義します。元のd0がd[0]と配列の添え字になります。

# ノード数
d = [784, 100, 50, 10]

#d0 = 784
#d1 = 100
#d2 = 50
#d3 = 10

階層数は、ノードの要素数-1になります。

# 階層数
layer = len(d) - 1

重み、バイアスには辞書を使いました。変数の定義を追加しています。

# 重み、バイアスの初期化
W = {}
b = {}
for i in range(layer):
    W[i+1] = he_normal(d[i],d[i+1])
    b[i+1] = np.zeros(d[i+1])

#W1 = he_normal(d0,d1)
#b1 = np.zeros(d1)
#W2 = he_normal(d1,d2)
#b2 = np.zeros(d2)
#W3 = he_normal(d2,d3)
#b3 = np.zeros(d3)

propagation関数のu,zにも辞書を使っています。0番目を入力層としています。

def propagation(layer, x, W, b):
    u = {}
    z = {}
    for i in range(layer+1):
        # 入力層
        if i == 0:
            z[i] = x
        # 中間層
        elif i < layer:
            u[i] = affine(z[i-1], W[i], b[i])
            z[i] = middle_func(u[i])
        # 出力層
        else:
            u[i] = affine(z[i-1], W[i], b[i])
            y = output_func(u[i])
    return u, z, y

    # 中間層(1)
    #u1 = affine(x, W1, b1)
    #z1 = middle_func(u1)
    # 中間層(2)
    #u2 = affine(z1, W2, b2)
    #z2 = middle_func(u2)
    # 出力層
    #u3 = affine(z2, W3, b3)
    #y = output_func(u3)
    #return u1, z1, u2, z2, u3, y

勾配計算は、関数にしました。同様に、du,dz,dW,dbにも辞書を利用しています。戻り値も辞書としています。

def back_propagation(layer, u, z, y, t, W, b):
    du = {}
    dz = {}
    dW = {}
    db = {}
    for i in range(layer, 0, -1):
        # 出力層
        if i == layer:
            du[i] = output_error_back(y, u[i], t)
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
        # 中間層
        else:
            du[i] = middle_back(dz[i], u[i], z[i])
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
    return du, dz, dW, db

        #du3 = output_error_back(y, u3, t[i:i+batch_size])
        #dz2, dW3, db3 = affine_back(du3, z2, W3, b3)
        #du2 = middle_back(dz2, u2, z2)
        #dz1, dW2, db2 = affine_back(du2, z1, W2, b2)
        #du1 = middle_back(dz1, u1, z1)
        #dz0, dW1, db1 = affine_back(du1, z0[i:i+batch_size], W1, b1)

重み、バイアスの調整部分の見直しです。

        # 重み、バイアスの調整
        for k in range(1, layer+1):
            W[k] = W[k] - eta*dW[k]
            b[k] = b[k] - eta*db[k]

        #W1 = W1 - eta*dW1
        #b1 = b1 - eta*db1
        #W2 = W2 - eta*dW2
        #b2 = b2 - eta*db2
        #W3 = W3 - eta*dW3
        #b3 = b3 - eta*db3

誤差の計算は、以下で行っています。

train_err = error_func(y_train, t_train)
test_err = error_func(y_test, t_test)

交差エントロピーのlogが0となる可能性があるので以下のように変更しました。

def cross_entropy_error(y, t):
    size = 1
    if y.ndim == 2:
        size = y.shape[0]
    #return -np.sum(t * np.log(y))/size
    return -np.sum(t * np.log(np.maximum(y, 1e-7)))/size

もう１点変更しています。データのシャッフルを事前に１回のみ実施していました。シャッフルはエポックごとに実施する必要があるようです。シャッフルをエポックのループの内側に移動します。

for j in range(epoch):
    # データのシャッフル(正解データも同期してシャフルする必要があるため一度、結合し分離)
    nx_t = np.concatenate([nx_train, t_train], axis=1)
    np.random.shuffle(nx_t)
    nx, t = np.split(nx_t, [nx_train.shape[1]], axis=1)

その他一部変数名の見直しを行いました。また、今後の対応を考慮し、勾配計算のmiddle_back関数に$ z $を渡すようにしました。
実行するためには、今までの関数も必要です。全体を一番下の（参考）に記載しましたので、そちらを先に実行してください。

import numpy as np
import time

# ノード数
d = [784, 100, 50, 10]
# 階層数
layer = len(d) - 1
# 重み、バイアスの初期化
W = {}
b = {}
for i in range(layer):
    W[i+1] = he_normal(d[i],d[i+1])
    b[i+1] = np.zeros(d[i+1])
# 学習率、バッチサイズ、エポックの設定
eta = 0.1
batch_size = 100
epoch = 50

# MNISTデータ読み込み
x_train, t_train, x_test, t_test = load_mnist('c:\\mnist\\')
# 入力データの正規化
nx_train = init_func(x_train)
nx_test  = init_func(x_test)

# 正解率、誤差初期化
train_rate = np.zeros(epoch+1)
test_rate = np.zeros(epoch+1)
train_err = np.zeros(epoch+1)
test_err = np.zeros(epoch+1)
# 実行(学習データ)
u_train, z_train, y_train = propagation(layer, nx_train, W, b)
train_rate[0] = accuracy_rate(y_train, t_train)
train_err[0] = error_func(y_train, t_train)
# 実行(テストデータ)
u_test, z_test, y_test = propagation(layer, nx_test, W, b)
test_rate[0] = accuracy_rate(y_test, t_test)
test_err[0] = error_func(y_test, t_test)
# 正解率、誤差表示
print(" 学習データ正解率 = " + str(train_rate[0]) + " テストデータ正解率 = " + str(test_rate[0]) +
      " 学習データ誤差 = " + str(train_err[0]) + " テストデータ誤差 = " + str(test_err[0]))

# 開始時刻設定
start_time = time.time()

for j in range(epoch):
    # データのシャッフル(正解データも同期してシャフルする必要があるため一度、結合し分離)
    nx_t = np.concatenate([nx_train, t_train], axis=1)
    np.random.shuffle(nx_t)
    nx, t = np.split(nx_t, [nx_train.shape[1]], axis=1)

    for i in range(0, nx.shape[0], batch_size):
        # 実行
        u, z, y = propagation(layer, nx[i:i+batch_size], W, b)
        # 勾配を計算
        du, dz, dW, db = back_propagation(layer, u, z, y, t[i:i+batch_size], W, b)
        
        # 重み、バイアスの調整
        for k in range(1, layer+1):
            W[k] = W[k] - eta*dW[k]
            b[k] = b[k] - eta*db[k]

    # 重み、バイアス調整後の実行(学習データ)
    u_train, z_train, y_train = propagation(layer, nx_train, W, b)
    train_rate[j+1] = accuracy_rate(y_train, t_train)
    train_err[j+1] = error_func(y_train, t_train)
    # 重み、バイアス調整後の実行(テストデータ)
    u_test, z_test, y_test = propagation(layer, nx_test, W, b)
    test_rate[j+1] = accuracy_rate(y_test, t_test)
    test_err[j+1] = error_func(y_test, t_test)
    # 正解率、誤差表示
    print(str(j+1) + " 学習データ正解率 = " + str(train_rate[j+1]) + " テストデータ正解率 = " + str(test_rate[j+1]) +
         " 学習データ誤差 = " + str(train_err[j+1]) + " テストデータ誤差 = " + str(test_err[j+1]))

# 終了時刻設定
end_time = time.time()
total_time = end_time - start_time
print("所要時間 = " +str(int(total_time/60))+" 分 "+str(int(total_time%60)) + " 秒")

実行条件

これからパラメータを変更しながら試していきます。以下の条件で実行します。

実行は、1回のみ
乱数は、seed値を設定し、毎回同じ乱数を発生するようにする。

np.random.seed(n)  # nは適当な数値にしてください。

乱数のseed値を変更すると、この結果通りとはならないかもしれませんが、十分参考にできると考えています。

まずは、基本のパラメータで実行してみました。
正解率は、以下のグラフのようになりました。実線がテストデータ正解率、破線が学習データ正解率です。

import numpy as np
import matplotlib.pyplot as plt
times = np.arange(0, epoch+1)
plt.plot(times, test_rate, label="Test Data", color="blue")
plt.plot(times, train_rate, label="Train Data", color="blue", linestyle="dashed")
plt.title("Accuracy rate")
plt.xlabel("epoch")
plt.ylabel("rate")
plt.legend()
plt.show()

学習データ誤差、テストデータ誤差についてもグラフ化します。

import numpy as np
import matplotlib.pyplot as plt
times = np.arange(0, epoch+1)
plt.plot(times, test_err, label="Test Error", color="blue")
plt.plot(times, train_err, label="Train Error", color="blue", linestyle="dashed")
plt.title("Error")
plt.xlabel("epoch")
plt.ylabel("value")
plt.legend()
plt.show()

また、パラメータ変更時の結果比較のため、以下のデータを表にします。

学習正解：50エポック実行後の学習データの正解率(%)
テスト正解：50エポック実行後のテストデータの正解率(%)
TS90%：テストデータの正解率が90%に達したエポック数
TS95%：テストデータの正解率が95%に達したエポック数
TS97%：テストデータの正解率が97%に達したエポック数
TR100%：学習データの正解率が100%に達したエポック数

学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
100.00	97.85	1	2	5	39

階層数・ノード数

階層数、ノード数を変更しながら試してみます。

１階層

まったく中間層がない場合です。活性化関数もありません。

# ノード数
d = [784, 10]

結果として、中間層がなくても、50エポックで、92.61%の正解率になりました。う～ん。中間層なくてもそれなりにできるじゃないですか、

２階層

中間層が１つです。中間層のノードが変わると結果がどのように変わるか見てみます。
1,5,10,50,100,500,1000で試してみます。

# ノード数
d = [784, 1, 10]

結果を表にします。所要時間、重みの数を追加します。
重みの数は、以下のように計算しています。
　ノード数100の場合
　　第１層： 784 $ \times $ 100 = 78400
　　第２層： 100 $ \times $ 10 = 1000
　　合計： 79400
参考として、1階層(中間層なし)も記載しました。

ノード数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間	重み数
1	37.55	37.67	-	-	-	-	1:28	794
5	90.81	89.94	21	-	-	-	1:46	3,970
参考	93.11	92.61	1	-	-	-	1:45	7,840
10	95.01	94.07	2	-	-	-	1:48	7,940
50	99.69	97.46	1	5	17	-	2:31	39,700
100	99.91	97.80	1	4	10	-	3:39	79,400
500	99.99	98.09	1	3	6	-	13:46	397,000
1000	99.99	98.20	1	3	7	-	24:31	794,000

1000ノードでは、私の非力なノートPC（Core i3、メモリ4GB)では、ほとんど固まってしまいました。

判明したこと

ノード数を増やすに従って、正解率も向上する。（あたりまえですが）
中間層のノード数が出力層のノード数より少ない場合は、１階層の場合より悪い結果となる。
ノード数が増えると急激にマシンリソースを消費する。

ノード数を増やせば増やすほど、重みの数が爆発的に増えます。精度が向上しますが、マシンリソースをかなり消費しました。階層数を増やすことで、リソースの消費を抑えて良い結果が出れば、その方がよいことになります。

３階層

中間層が２つになります。中間層はどのような組み合わせがよいか見てみます。各層10,50,100,500で試してみます。

# ノード数
d = [784, 10, 10, 10]

ノード数1	ノード数2	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間	重み数
10	10	95.44	93.98	1	-	-	-	1:44	8,040
10	50	97.13	95.58	1	9	-	-	1:54	8,840
10	100	97.81	95.61	1	11	-	-	2:21	9,840
10	500	97.33	95.32	2	17	-	-	4:23	17,840
50	10	100.00	97.39	1	3	9	-	2:49	39,800
50	50	100.00	97.74	1	3	8	-	2:58	42,200
50	100	100.00	97.72	1	2	7	48	3:27	45,200
50	500	100.00	97.63	1	2	7	39	5:45	69,200
100	10	100.00	97.85	1	3	6	-	3:48	79,500
100	50	100.00	97.85	1	2	5	39	3:55	83,900
100	100	100.00	97.99	1	2	6	38	3:57	89,400
100	500	100.00	97.92	1	2	5	33	6:12	133,400
500	10	100.00	98.19	1	2	5	-	13:51	397,100
500	50	100.00	98.27	1	2	4	25	13:45	417,500
500	100	100.00	98.30	1	2	5	31	15:47	443,000
500	500	100.00	98.17	1	2	4	27	25:17	647,000

２層の1000ノードの場合の正解率を上回ったのは、500-50,500-100です。
２層にした効果は出ているようです。正解率が97%に達するエポック数はかなり早くなっていますが、最終的な正解率の向上は思ったほどよくなっていません。

判明したこと

第１層のノード数が小さいと第２層を大きくしても正解率は上がらない。
第１層のノード数への依存度が大きく、第２層のノード数を増やしてもあまり正解率が変化しない。

多階層

今度は、ノード数を一定にして階層数を深くしてみます。

２階層で各層のノード数が100の場合の設定例

# ノード数
d = [784, 100, 10]

100ノードの場合の結果です。参考として、２階層、1000ノードの場合を追加しています。

階層数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間	重み数
2	99.91	97.80	1	4	10	-	3:39	79,400
3	100.00	97.99	1	2	6	38	4:10	89,400
4	100.00	97.99	1	2	5	29	4:47	99,400
5	100.00	98.03	1	2	5	31	5:04	109,400
6	100.00	98.07	1	1	4	29	6:06	119,400
7	100.00	97.97	1	1	5	34	6:50	129,400
(参考)	99.99	98.20	1	3	7	-	24:31	794,000

100ノードで階層数を単純に増やしただけでは、2層の1000ノードにはおよびませんでした。

判明したこと

階層数を増やすと、テストデータ正解率は向上したが、向上率はわずかであった。
階層数を増やすと、一定の正解率には早く達した。

単に、階層数を増やせばよいわけではないようです。

今回、誤差も表示するように改造しましたので、学習データの誤差、テストデータの誤差を見て見ます。参考として、２階層、1000ノードの場合を追加しています。

階層数	学習正解	テスト正解	学習誤差	テスト誤差
2	99.91	97.80	0.0107	0.0770
3	100.00	97.99	0.0011	0.0919
4	100.00	97.99	0.0003	0.1091
5	100.00	98.03	0.0001	0.1083
6	100.00	98.07	0.0001	0.1157
7	100.00	97.97	0.0000	0.1401
(参考)	99.99	98.20	0.0053	0.0589

勾配は、学習データの誤差を小さくするよう計算するのでした。学習データに対する誤差が非常に小さいため学習が進まなくなっていると想像します。

2階層1000ノードの場合の重みの数は、79万個です。それに対して、3階層100ノード、100ノードの場合の重みの数は、約10分の1の8.9万個です。ノード数を増やせば正解率はあがりましたが、それだけマシンリソースを食います。層を増やした方がリソースを食わないことがわかりました。

重みの初期値

重みの初期値を変更して試してみます。

初期化関数

重みの初期値を、he_normalから、he_uniform、lecun_normal、lecun_uniform、glorot_normal、glorot_uniformに変更してみます。

for i in range(layer):
    W[i+1] = he_normal(d[i],d[i+1])
    b[i+1] = np.zeros(d[i+1])

初期化関数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
he_normal	100.00	97.85	1	2	5	39
he_uniform	100.00	97.85	1	2	6	36
lecun_normal	100.00	98.07	1	3	6	38
lecun_uniform	100.00	98.06	1	2	6	41
glorot_normal	100.00	97.91	1	2	5	37
glorot_uniform	100.00	97.99	1	2	7	41

lecunが一番よい結果となりました。乱数のため、seed値を変更すると別の結果になるかもしれません。

それぞれの重み、バイアスの値をヒストグラムにしてみます。
初期値のヒストグラム表示です。バイアスは、0のため省略しています。

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(10,2.2))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
# 重みのヒストグラム
for i in range(1, layer+1):
    fW = W[i].flatten()
    plt.subplot(1, 3, i)
    plt.title("W[" + str(i) + "]")
    plt.hist(fW,bins=50) 

plt.show()

50エポック後のヒストグラム表示です。

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(10,5))
plt.subplots_adjust(wspace=0.3, hspace=0.3)
# 重みのヒストグラム
for i in range(1, layer+1):
    fW = W[i].flatten()
    plt.subplot(2, 3, i)
    plt.title("W[" + str(i) + "]")
    plt.hist(fW,bins=50)
# バイアスのヒストグラム
for i in range(1, layer+1):
    fb = b[i].flatten()
    plt.subplot(2, 3, i+3)
    plt.title("b[" + str(i) + "]")
    plt.hist(fb,bins=50)  

plt.show()

以下のようなグラフとなりました。正解率はどれもほぼ同じですが、50エポック後の値の分布はばらばらです。同じ値に収束するわけではないことがわかりました。よく見ると、初期値と似た傾向があります。初期値によって、収束する値が変わっていくことがわかりました。

he_normal

初期値

50エポック後

he_uniform

初期値

50エポック後

lecun_normal

初期値

50エポック後

lecun_uniform

初期値

50エポック後

glorot_normal

初期値

50エポック後

glorot_uniform

初期値

50エポック後

固定値

次に、オール0とオール1で試し見ます。

# オール0
for i in range(layer):
    W[i+1] = np.zeros((d[i],d[i+1]))
    b[i+1] = np.zeros(d[i+1])
# オール1
for i in range(layer):
    W[i+1] = np.ones((d[i],d[i+1]))
    b[i+1] = np.zeros(d[i+1])

初期化関数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
np.zeros	11.24	11.35	-	-	-	-
np.ones	11.24	11.35	-	-	-	-

オール0、オール1とも収束しませんでした。
$ u_i $は、以下のように求めました。

u_i = w_{i1}z_1 + w_{i2}z_2 + \cdots + w_{in}z_n + b_i

オール0の場合、重み、バイアスとも0のため$ u_i $は、すべて0になります。重みの勾配は、$ z_i $を掛けることで計算しました。ReLU関数を利用しているため$ z_i $も0になります。よって、重みは、0のまま学習はされません。
オール1の場合、$ u_i $は、$ z_i $をすべて足した値になります。$ u_i $は、すべて同じ値です。$ z_i $もすべて同じ値のため、重みの勾配もすべて同じ値になります。よって、各層の重みがすべて同じ値になり、学習されません。

乱数

初期値は、固定値だとうまくいきませんでした。では、乱数であれば何でもよいか確認してみます。
まずは、値の範囲を指定して乱数を生成する関数を作成します。

def random_w(d_1, d, min=0, max=1):
    return np.random.rand(d_1,d) *  (max-min) + min

for i in range(layer):
    # 0～1の乱数
    W[i+1] = random_w(d[i],d[i+1])
    # -1～1の乱数
    W[i+1] = random_w(d[i],d[i+1],-1,1)
    # -0.1～0.1の乱数
    W[i+1] = random_w(d[i],d[i+1],-0.1,0.1)

0～1、-1～1、-0.1～0.1の乱数で以下の結果となりました。

乱数範囲	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
0～1	11.24	11.35	-	-	-	-
-1～1	97.48	95.28	6	34	-	-
-0.1～0.1	100.00	97.87	1	4	9	45

-0.1～0.1の場合のヒストグラムを表示します。

初期値

50エポック後

0～1の場合は、固定値と同じでまったく収束しませんでした。重みがすべてプラスの値のため、結果として勾配計算がうまくいかないのか、負の値まで拡張すると、それなりの結果となりました。

判明したこと

重みの初期値を定数にすると学習が進まない。
乱数を利用すると学習が進む。
特に、Heの初期値などのようにノード数により重みの初期値を調整すると良い。

重みの初期値は、学習に際して、非常に重要だということがわかりました。

バイアス

バイアスも重みと同様に固定値、乱数を初期値として試してみました。
重みと同様に、値の範囲を指定して乱数を生成する関数を作成します。

def random_b(d, min=0, max=1):
    return np.random.rand(d) *  (max-min) + min

for i in range(layer):
    # 0(固定)
    b[i+1] = np.zeros(d[i+1])
    # 1(固定) 
    b[i+1] = np.ones(d[i+1])
    # 0～1(乱数)
    b[i+1] = random_b(d[i+1])
    # -1～1(乱数)
    b[i+1] = random_b(d[i+1],-1,1)
    # -0.1～0.1(乱数)
    b[i+1] = random_b(d[i+1],-0.1,0.1)

初期値	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
0(固定)	100.00	97.85	1	2	5	39
1(固定)	100.00	97.87	1	3	9	43
0～1(乱数)	100.00	98.05	1	2	6	37
-1～1(乱数)	100.00	97.90	1	2	5	48
-0.1～0.1(乱数)	100.00	98.02	1	2	5	37

0～1(乱数)の場合のバイアスの変化をヒストグラムで見てみましょう。

初期値

50エポック後

学習後のバイアスもほぼ0～1に収まっています。比較のため、-0.1～0.1(乱数)のヒストグラムも見てみましょう。

初期値

50エポック後

バイアスも初期値に依存し、学習後の値が変わることがわかりました。

判明したこと

バイアスの初期値を変えても、重みに比べ正解率にあまり変化がない。

正規化関数

今までは、学習データ、テストデータ個別にmin_max(0～1)、Zスコア(平均0、分散1)を計算していました。本来であればテストデータにも学習データで求めた、最小値・最大値(min_max)、平均・分散(z_score)を利用すべきでした。テストデータの情報は利用してはいけないのでした。
MNISTでは、学習データ、テストデータとも最小値=0、最大値=255のため、今までのmin_maxの結果に影響はありませんでした。よかった。
実装の変更をします。学習時、最小値・最大値(min_max)、平均・分散(z_score)を計算し、テスト時は、学習データで求めた値を利用できるようにします。学習時は、パラメータを指定せずに、計算を実施、テスト時には、学習で求めた値をパラメータして渡し利用します。

# 正規化関数
def min_max(x, x_min=None, x_max=None, axis=None):
    if x_min is None:
        x_min = np.min(x, axis=axis, keepdims=True) # 最小値を求める
    if x_max is None:
        x_max = np.max(x, axis=axis, keepdims=True) # 最大値を求める
    return (x-x_min)/np.maximum((x_max-x_min),1e-7), x_min, x_max

def z_score(x, x_mean=None, x_std=None, axis=None):
    if x_mean is None:
        x_mean = np.mean(x, axis=axis, keepdims=True) # 平均値を求める
    if x_std is None:
        x_std  = np.std(x, axis=axis, keepdims=True)  # 標準偏差を求める
    return (x-x_mean)/np.maximum(x_std,1e-7), x_mean, x_std

min_maxの場合は、以下となります。

# 入力データの正規化
nx_train, train_min, train_max = min_max(x_train)
nx_test,  test_min,  test_max  = min_max(x_test, train_min, train_max)

z_scoreの場合は、以下となります。

# 入力データの正規化
nx_train, train_mean, train_std = z_score(x_train)
nx_test,  test_mean,  test_std  = z_score(x_test, train_mean, train_std)

また、本来であれば、データ系列ごとに正規化する必要があります。今までは、全データを対象にしていました。MNISTでは、すべてのデータが0～255のピクセル値のためそれでも問題ありませんでした。
今回は、軸を指定して実行してみます。

データ系列ごとに正規化するためには、axis=0をパラメータとして指定します。（z_scoreの例）

# 入力データの正規化
nx_train, train_mean, train_std = z_score(x_train, axis=0)
nx_test,  test_mean,  test_std  = z_score(x_test, train_mean, train_std, axis=0)

全データ(axis指定なし）、データ系列ごと(axis=0)で実行してみます。

関数	axis	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
mix_max	なし	100.00	97.85	1	2	5	39
mix_max	0	100.00	97.83	1	2	5	40
z_score	なし	100.00	97.96	1	1	4	23
z_score	0	100.00	97.41	1	2	6	32
なし	-	11.24	11.35	-	-	-	-

mix_max(axis=0)、z_score(axis=0)の場合の12エポックまでをグラフ化してみます。

序盤は、z_scoreの方が、学習が早く進むことがわかりました。
Zスコア(なし)は、なんと、1エポックで正解率が95%を超えました。今までで最高です。また、学習データの正解率が100%に達するまでのエポック数23も今までで最高です。学習を早める効果は、かなりありそうです。

判明したこと

データを正規化することで、学習が早く進む。

活性化関数(中間層)

活性化関数を変更する場合、変更した活性化関数の勾配を求める関数を求める必要があります。前回の法則を利用し勾配関数を決めながら確認していきましょう。
前回の法則です。

**　法則１：勾配は、ネットワークの各線の勾配の積となる。**
**　法則２：分岐の場合は、それぞれの勾配の和となる。**
**　法則３：複雑な関数は分解し、法則１、法則２を利用し勾配を求める。**
**　（$ + $定数）：定数を加えた場合の勾配 $ =1 $**
**　（$ \times $定数）：定数を掛けた場合の勾配 $ = $ 定数**
**　（$ \sum $）：$ \sum $の勾配 $ =1 $**
**　（$ \times $）：２つの変数の掛け算の勾配は、それぞれ逆の変数となる**
**　（関数）：関数の勾配は、数学の力を借りる**
　　　　　　$ \log(x) $ → $ \frac{1}{x} $、$ \exp(x) $ → $ \exp(x) $、$ \frac{1}{x} $ → $ -\frac{1}{x^2} $

定義済み活性化関数

sigmoid

sigmoidは、以下の関数でした。

f(x) = \frac{1}{1+\exp(-x)}

関数を分解すると、以下です。
第１層： $ -1 $を掛ける。
　　　$ x $ → $ -x $
第２層： $ \exp $
　　　$ -x $ → $ \exp(-x) $
第３層： $ +1 $する。
　　　$ \exp(-x) $ → $ 1 + \exp(-x) $
第４層：逆数をとる。
　　　$ 1 + \exp(-x) $ → $ \frac{1}{1 + \exp(-x)} $
ネットワーク図にしてみます。

勾配の計算は、今回は、前から求めます。
第１層： $ -1 $を掛ける。
　　　法則３（$ \times $定数）より、掛け算の勾配は、掛けた定数の$ -1 $
第２層： $ \exp $
　　　法則３（関数）より、勾配は、そのまま、$ \exp(-x) $　
第３層： $ +1 $する。
　　　法則３（$ + $定数）より、足し算の勾配は、$ 1 $固定
第４層：逆数をとる。
　　　$ 1 + \exp(-x) $　の逆数のため、法則３（関）より、マイナス２乗分の１で、勾配は、$ -\frac{1}{(1 + \exp(-x))^2} $　
これらをすべて掛けると、以下になります。

\begin{align}
f^{'}(x) &= -1 \times \exp(-x) \times 1 \times -\frac{1}{(1 + \exp(-x))^2}\\
&= \frac{\exp(-x)}{(1 + \exp(-x))^2}
\end{align}

このままでもよいのですが、ここで技を使うようです。分子に $ +1 $、$ -1 $します。足して、引けば同じですものね。そして、分子を$ 1 + \exp(-x) $と$ -1 $に、ばらして書いてみます。

\begin{align}
f^{'}(x) &= \frac{\exp(-x)}{(1 + \exp(-x))^2}\\
&= \frac{1 + \exp(-x) - 1}{(1 + \exp(-x))^2}\\
&= \frac{1 + \exp(-x)}{(1 + \exp(-x))^2}-\frac{1}{(1 + \exp(-x))^2}\\
&= \frac{1}{1 + \exp(-x)}-\frac{1}{(1 + \exp(-x))^2}
\end{align}

よく見ると$ f(x) $が含まれています。

f(x) = \frac{1}{1+\exp(-x)}

よって

\begin{align}
f^{'}(x) &= \frac{1}{1 + \exp(-x)}-\frac{1}{(1 + \exp(-x))^2}\\
&= f(x) - f(x)^2
\end{align}

勾配関数を実装してみます。

def sigmoid_back(z):
    return z - z**2

活性化関数の勾配のパラメータは、$ x $でした。sigmoidの勾配は、$ z $で計算するため、middle_backに$ z $をパラメータとして追加が必要です。先に、パラメータを追加したのはこのためです。

勾配関数は、次のように定義しました。

def middle_func(x):
    return sigmoid(x)
def middle_back_func(x, z):
    return sigmoid_back(z)

tanh

勾配の計算はあきらめて、数学の力をかりましょう。

\begin{align}
f(x) &= \tanh(x)\\
f^{'}(x) &= \frac{1}{\cosh(x)^2}
\end{align}

def tanh_back(x):
    return 1/np.cosh(x)**2

def middle_func(x):
    return tanh(x)
def middle_back_func(x, z):
    return tanh_back(x)

identity

恒等関数です。活性化関数を通さない場合です。勾配は説明するまでもなく、$ 1 $ですね。

def identity_back(x):
    return 1

def middle_func(x):
    return identity(x)
def middle_back_func(x, z):
    return identity_back(x)

softplus

他に活性化関数として利用される関数でも試してみました。
まずは、softplus関数です。ReLUを少し滑らかにしたグラフになります。青がsoftplus、オレンジがReLUになります。

以下の式で表します。

f(x) = \log(1 + \exp(x))

以下のネットワーク図になります。

勾配は、以下の計算になります。

\begin{align}
f^{'}(x) &= \exp(x) \times 1 \times \frac{1}{1 + \exp(x)}\\
         &= \frac{\exp(x)}{1 + \exp(x)}
\end{align}

ここで、$ e^{-x} = \frac{1}{e^x} $、$ e^x \times e^y = e^{(x+y)} $の性質を利用して、変形します。

\begin{align}
f^{'}(x) &= \frac{\exp(x)}{1 + \exp(x)}\\
         &= \frac{1}{\exp(-x)}\frac{1}{(1 + \exp(x))}\\
         &= \frac{1}{\exp(-x) + \exp(-x)\exp(x)}\\
         &= \frac{1}{\exp(-x) + \exp(0)}\\
         &= \frac{1}{1 + \exp(-x)}
\end{align}

どこかで見た式になりました。そう、sigmoid関数です。softplusの勾配は、sidmoidになります。

def softplus(x):
    return np.log(1+np.exp(x))
def softplus_back(x):
    return 1/(1 + np.exp(-x))

def middle_func(x):
    return softplus(x)
def middle_back_func(x, z):
    return softplus_back(x)

softsign

次は、softsignです。tanh関数を滑らかにしたグラフになります。青がsoftsign、オレンジがtanhになります。

以下の式で表します。

f(x) = \frac{x}{1 + |x|}

勾配は、0以上の場合と0未満の場合に分けて考えてみましょう。

f(x) = \left\{
\begin{array}{ll}
\frac{x}{1 + x} & (x \geq 0)\\
\frac{x}{1 - x} & (x \lt 0)
\end{array}
\right.

勾配のネットワーク図、勾配の計算は、それぞれ以下のようになります。

$ x \geq 0 $の場合、$ |x| = x $

\begin{align}
f^{'}(x) &= \frac{1}{1 + x} + 1 \times -\frac{1}{(1 + x)^2} \times x\\
         &= \frac{1}{1 + x} - \frac{x}{(1 + x)^2}\\
         &= \frac{1+x}{(1 + x)^2} - \frac{x}{(1 + x)^2}\\
         &= \frac{1}{(1 + x)^2}
\end{align}

$ x \lt 0 $の場合、$ |x| = -x $

\begin{align}
f^{'}(x) &= \frac{1}{1 - x} + (-1) \times 1 \times -\frac{1}{(1 - x)^2} \times x\\
         &= \frac{1}{1 - x} + \frac{x}{(1 - x)^2}\\
         &= \frac{1-x}{(1 - x)^2} + \frac{x}{(1 - x)^2}\\
         &= \frac{1}{(1 - x)^2}
\end{align}

勾配計算の説明を省略していますが、大丈夫でしょうか？掛け算は、お互いの逆の値になるのでした。逆数の計算は、$ -\frac{1}{x^2} $でした。
最終的に勾配は、両方とも以下になります。

f^{'}(x) = \frac{1}{(1 + |x|)^2}

def softsign(x):
    return x/(1+np.absolute(x))
def softsign_back(x):
    return 1/(1+np.absolute(x))**2

def middle_func(x):
    return softsign(x)
def middle_back_func(x, z):
    return softsign_back(x)

step

最後にステップ関数です。ReLUの勾配関数です。
勾配は、$ x=0 $の場合は、正確には定まりませんが、$ 0 $ですよね。

\begin{align}
f(x) &= \left\{
\begin{array}{ll}
1 & (x \gt 0)\\
0 & (x \leq 0)
\end{array}
\right.
\\
f^{'}(x) &= 0
\end{align}

def step(x):
    return np.where(x > 0, 1, 0)
def step_back(x):
    return x*0

def middle_func(x):
    return step(x)
def middle_back_func(x, z):
    return step_back(x)

結果

tanhが一番良い結果となりました。sigmoidはあまりよくありません。恒等関数、すなわち、活性化関数なしが以外に健闘しています。

活性化関数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
relu	100.00	97.85	1	2	5	39
sigmoid	98.35	97.38	4	18	42	-
tanh	100.00	97.89	1	3	7	45
identity	93.31	92.27	1	-	-	-
softplus	99.87	97.68	1	5	15	-
softsign	99.97	97.79	1	4	11	-
step	52.70	52.40	-	-	-	-

正解率をグラフにしてみます。実線がテストデータ正解率、破線が学習データ正解率です。

ReLU関連

ReLUの負の部分(0固定部分)に一定の傾きを持った関数も利用されるようです。
以下の式になります。

f(x) =  \left\{
\begin{array}{ll}
x & (x \gt 0) \\
\alpha x & (x \leq 0)
\end{array}
\right.

$ \alpha $は、$ 0 \lt \alpha \lt 1 $です。$ \alpha = 0 $の場合はReLU、$ \alpha = 1 $の場合は恒等関数です。
$ \alpha = 0.1 $のグラフは、以下のようになります。

以下の方法があるようです。

Leaky ReLU
　　$ \alpha $を定数とする。
RReLU
　　$ \alpha $をランダムに与える。
PReLU
　　$ \alpha $を重みと同じようにパラメータとらえ、学習の対象とする。

Leaky ReLU

Leaky ReLUを実装します。関数は、以下になります。勾配は、負の場合、$ \alpha $ですよね。

def leaky_relu(x, alpha):
    return np.maximum(x * alpha, x)

def leaky_relu_back(x, alpha):
    return np.where(x > 0, 1, alpha)

$ \alpha $を変更しながら確認するため、活性化関数、勾配関数にパラメータを渡せるように実装を変更します。propagation、back_propagationにパラメータを追加します。また、propagation、back_propagationを呼び出している部分にもパラメータを追加します。

def propagation(layer, x, W, b, param):
    u = {}
    z = {}
    for i in range(layer+1):
        # 入力層
        if i == 0:
            z[i] = x
        # 中間層
        elif i < layer:
            u[i] = affine(z[i-1], W[i], b[i])
            z[i] = middle_func(u[i], param)
        # 出力層
        else:
            u[i] = affine(z[i-1], W[i], b[i])
            y = output_func(u[i])
    return u, z, y

def back_propagation(layer, u, z, y, t, W, b, param):
    du = {}
    dz = {}
    dW = {}
    db = {}
    for i in range(layer, 0, -1):
        # 出力層
        if i == layer:
            du[i] = output_error_back(y, u[i], t)
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
        # 中間層
        else:
            du[i] = middle_back(dz[i], u[i], z[i], param)
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
    return du, dz, dW, db

def middle_back(dz, x, z, param):
    du = middle_back_func(x, z, param)
    return dz * du

パラメータは、汎用的に利用できるように辞書を利用しました。パラメータ設定および受け渡し部分は、以下のようにしました。

def middle_func(x, param):
    return leaky_relu(x, param["alpha"])
def middle_back_func(x, z, param):
    return leaky_relu_back(x, param["alpha"])

# パラメータ設定
param = {}
param["alpha"] = 0.1

結果

パラメータを0.1～0.9まで0.1ごとに実行してみました。参考のため、0.0=ReLU、1.0=恒等関数も含めます。

$ \alpha $	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
0.0	100.00	97.85	1	2	5	39
0.1	100.00	97.91	1	2	6	45
0.2	100.00	97.94	1	3	7	45
0.3	100.00	97.92	1	3	7	-
0.4	99.97	97.65	1	4	10	-
0.5	99.55	97.36	1	4	14	-
0.6	99.18	97.34	1	5	23	-
0.7	98.62	97.02	1	8	32	-
0.8	97.00	95.48	1	20	-	-
0.9	94.83	93.76	1	-	-	-
1.0	93.31	92.27	1	-	-	-

正解率をグラフにしてみます。実線がテストデータ正解率、破線が学習データ正解率です。

$ \alpha $が0.3までは、ReLUを上回りました。ただし、大きな差はないようです。$ \alpha $が0.5より大きいと良い結果は得られませんでした。

RReLU、PReLUは、実装できるレベルまで理解できていないため、今回は、試していません。
また、maxout関数が良いとの情報も見ましたが、こちらも実装できるレベルまで理解が及んでいないため見送ります。

その他関数

せっかくなので思いつく関数で試してみます。

２乗

説明は、不要だと思いますが、勾配は、$ 2x $です。

\begin{align}
f(x) &= x^2\\
f^{'}(x) &= 2x
\end{align}

def power2(x):
    return x**2
def power2_back(x):
    return 2 * x

def middle_func(x):
    return power2(x)
def middle_back_func(x, z):
    return power2_back(x)

オーバフローが発生したため学習率を0.01に変更し実行しました。

eta = 0.01

３乗

勾配は、$ 3x^2 $です。

\begin{align}
f(x) &= x^3\\
f^{'}(x) &= 3x^2
\end{align}

def power3(x):
    return x**3
def power3_back(x):
    return 3 * x**2

def middle_func(x):
    return power3(x)
def middle_back_func(x, z):
    return power3_back(x)

オーバフローが発生したため学習率を0.001に変更し実行しました。

eta = 0.001

平方根

次にルートです。0以下は、$ 0 $とします。

f(x) = \left\{
\begin{array}{ll}
\sqrt{x} & (x \gt 0)\\
0 & (x \leq 0)
\end{array}
\right.

ルートの勾配は、数学の力を借りて、$ \frac{1}{2 \sqrt{x}}　$とします。

f^{'}(x) = \left\{
\begin{array}{ll}
\frac{1}{2 \sqrt{x}} & (x \gt 0)\\
0 & (x \leq 0)
\end{array}
\right.

def sqrt(x):
    return np.where(x > 0, np.sqrt(np.maximum(x, 1e-7)), 0)
def sqrt_back(x):
    return np.where(x > 0, 1/(2 * np.sqrt(np.maximum(x, 1e-7))), 0)

def middle_func(x):
    return sqrt(x)
def middle_back_func(x, z):
    return sqrt_back(x)

2乗(0以下は0)

2乗の成績が予想以上だったので、ReLUのように0以下は、0として試してみます。

f(x) = \left\{
\begin{array}{ll}
x^2 & (x \gt 0)\\
0 & (x \leq 0)
\end{array}
\right.

f^{'}(x) = \left\{
\begin{array}{ll}
2x & (x \gt 0)\\
0 & (x \leq 0)
\end{array}
\right.

def power2_0(x):
    return np.where(x > 0, x**2, 0)
def power2_0_back(x):
    return np.where(x > 0, 2 * x, 0)

def middle_func(x):
    return power2_0(x)
def middle_back_func(x, z):
    return power2_0_back(x)

オーバフローが発生したため学習率を0.01に変更し実行しました。

eta = 0.01

絶対値

f(x) = |x|

勾配は、$ x \gt 0 $の場合、$ 1 $、$ x \lt 0 $の場合、$ -1 $ですね。$ 0 $の場合は、不定ですが、$ 0 $とします。

f^{'}(x) = \left\{
\begin{array}{ll}
1 & (x \gt 0)\\
0 & (x = 0 )\\
-1 & (x \lt 0)
\end{array}
\right.

def absolute(x):
    return np.absolute(x)
def absolute_back(x):
    return np.where( x > 0, 1, np.where( x < 0, -1, 0))

def middle_func(x):
    return absolute(x)
def middle_back_func(x, z):
    return absolute_back(x)

sin

意味があるかどうかわかりませんが、sin,cosも確認してみます。勾配は、数学の力を借ります。

\begin{align}
f(x) &= \sin(x)\\
f^{'}(x) &= \cos(x)
\end{align}

def sin(x):
    return np.sin(x)
def sin_back(x):
    return np.cos(x)

def middle_func(x):
    return sin(x)
def middle_back_func(x, z):
    return sin_back(x)

学習率を0.01に変更し実行しました。

eta = 0.01

cos

\begin{align}
f(x) &= \cos(x)\\
f^{'}(x) &= -\sin(x)
\end{align}

def cos(x):
    return np.cos(x)
def cos_back(x):
    return -np.sin(x)

def middle_func(x):
    return cos(x)
def middle_back_func(x, z):
    return cos_back(x)

学習率を0.01に変更し実行しました。

eta = 0.01

結果

どうも勾配が大きくなるとオーバフローするようです。学習率を小さくすることで回避しました。そう言えば、通常利用される活性化関数の勾配は、0～1の間でした。
どんな関数でも学習率を選べばそれなりの結果が得られました。

活性化関数	学習率	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
２乗	0.01	99.98	97.49	1	4	13	-
３乗	0.001	96.77	94.49	9	-	-	-
平方根	0.1	98.13	96.65	1	6	-	-
２乗(0以下は0)	0.01	99.72	97.15	1	6	30	-
絶対値	0.1	100.00	97.82	1	1	25	22
sin	0.01	100.00	96.63	3	8	-	-
cos	0.01	100.00	96.80	3	8	-	-

活性化関数（出力層）＋損失関数

分類問題は、出力層の活性化関数はソフトマックス関数、損失関数は交差エントロピー関数が一般的なようです。ここでは、回帰問題で利用される二乗和誤差関数を利用し確認してみましょう。

恒等関数＋二乗和誤差

二乗和誤差の勾配を求めます。二乗和誤差は、以下の関数でした。

E = \frac{1}{2}\sum_k(y_k-t_k)^2

ネットワーク図を描いてみます。

$ y_1 $の勾配を計算します。

1 \times 2(y_1 - t_1) \times 1 \times \frac{1}{2} = y_1 - t_1

一般的には、$ y_i - t_i $になります。どこかで見た式ですね。そうです。ソフトマックス＋交差エントロピーの勾配と同じです。
恒等関数の勾配は、$ 1 $なので、出力層の活性化関数＋損失関数の勾配は、$ y_i - t_i $です。

def identity_mean_squared_error_back(y, x, t):
    size = 1
    if y.ndim == 2:
        size = y.shape[0]
    return (y - t)/size

def output_func(x):
    return identity(x)

def error_func(y, t):
    return mean_squared_error(y, t)

def output_error_back_func(y, x, t):
    return identity_mean_squared_error_back(y, x, t)

ソフトマックス＋二乗和誤差

ソフトマックスと二乗和誤差の組み合わせにしてみます。
ソフトマックスの勾配は、$ f^{'}(y_i)y_i - (\sum_kf^{'}(y_k)y_k)y_i $でした。ここに、二乗和誤差の勾配、$ y_i - t_i $を代入します。

(y_i - t_i)y_i - (\sum_k(y_k - t_k)y_k)y_i

実装してみます。

def softmax_mean_squared_error_back(y, x, t):
    if y.ndim == 1:
        size = 1
        sum_y = np.sum((y-t)*y)
    elif y.ndim == 2:
        size = y.shape[0]
        sum_y = np.sum((y-t)*y, axis=1, keepdims=True)
    return ((y-t)*y-sum_y*y)/size

def output_func(x):
    return softmax(x)

def error_func(y, t):
    return mean_squared_error(y, t)

def output_error_back_func(y, x, t):
    return softmax_mean_squared_error_back(y, x, t)

結果

これで何が言えるかわかりませんが、分類問題のため、出力層の活性化関数は、ソフトマックス、損失関数は交差エントロピーを利用します。

出力層関数	損失関数	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
ソフトマックス	交差エントロピー	100.00	97.85	1	2	5	39
恒等関数	二乗和誤差	98.98	97.55	1	5	19	-
ソフトマックス	二乗和誤差	99.30	97.51	2	8	21	-

学習率

学習率を 10, 1、0.1、0.01、0.001、0.0001に変更し試してみました。

結果

学習率が学習に大きな影響があることがわかりました。

学習率	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%
10.0	11.23	11.35	-	-	-	-
1.0	98.56	95.92	1	4	-	-
0.1	100.00	97.85	1	2	5	39
0.01	98.08	97.08	3	15	44	-
0.001	92.05	92.37	21	-	-	-
0.0001	78.69	79.22	-	-	-	-

学習率1.0のグラフを見てみましょう。

正解率が上下しなかなか収束しません。
ここで、各層の重みから３個づつ取り出してエポックごとの変化を見てみましょう。

重みの値が上下し収束に向かっていません。

次に、学習率0.001のグラフを見てみましょう。

今度は、学習率が徐々にしか向上しません。
重みの変化も見てみましょう。

学習率が1.0の場合と傾向が全く違いますね。

判明したこと

学習率が大きければ、重みが収束しない。
学習率が小さければ、正解率がなかなか収束しない。

バッチサイズ

バッチサイズを1(オンライン学習)、10、50、100、500、1000、60000(バッチ学習)に変更して試してみます。

結果

バッチサイズが50の場合が最高でした。
バッチサイズが1の場合、収束しない上、時間も非常にかかりました。

バッチサイズ	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間	バッチ数
1	11.24	11.35	-	-	-	-	40:08	3000000
10	99.50	97.34	1	1	4	-	7:52	300000
50	100.00	97.90	1	2	5	26	3:57	60000
100	100.00	97.85	1	2	5	39	3:55	30000
500	99.31	97.53	2	7	23	-	3:47	6000
1000	98.10	97.15	3	15	43	-	3:34	3000
60000	84.24	84.92	-	-	-	-	4:30	50

バッチサイズが1の場合、学習率を0.001に変更し実行すると良い結果を得られました。学習率が大きいとその画像にフィットするように重みが調整されるためと考えられます。バッチサイズに合わせて学習率を変更する必要がありそうです。

バッチサイズ	学習率	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間	バッチ数
1	0.001	100.00	97.86	1	2	5	37	43:52	3000000
10	0.01	100.00	97.88	1	2	5	37	7:18	300000

判明したこと

バッチサイズに合わせて学習率を調整する必要がある。
バッチサイズが小さすぎると時間がかかる。
バッチサイズが大きいとなかなか収束しない。

その他、試したこと

重みを平均の画像の値に

重みの初期値を各数字の平均値とした場合、どのような値になるか試してみました。
まずは、各数字の平均の画像の値を求めます。平均の画像を表示してみます。

import matplotlib.pyplot as plt
import numpy as np
# MNISTデータ読み込み
x_train, t_train, x_test, t_test = load_mnist('c:\\mnist\\')
# 各数字の個数
num_cnt = np.sum(t_train,axis=0)
for i in range(10):
    # 各数字のデータのみに変換
    num_data = (x_train.T*t_train[:,i].T).T
    # 平均のデータを計算
    data = np.sum(num_data,axis=0)/num_cnt[i]
    # 画像表示
    plt.subplot(2, 5, i+1)
    plt.imshow(data.reshape(28,28), 'gray')
plt.show()

おお、ちゃんと各数字に見えますね。
次に、768×10の中間層なしの場合を考えます。重みは、768×10個になります。重みを各数字の平均値に設定します。

# ノード数
d = [784, 10]

# 重みの設定
# 各数字の個数
num_cnt = np.sum(t_train,axis=0)
for i in range(10):
    # 各数字のデータのみに変換
    num_data = (x_train.T*t_train[:,i].T).T
    # 平均のデータを計算
    data = np.sum(num_data,axis=0)/num_cnt[i]
    #　重みを設定
    W[1][:,i] = data

学習なしで、以下の結果になりました。
　学習データ正解率　：62.33%
　テストデータ正解率：63.10%
もう少し良い値になるかと思っていたのですが、学習なしでもある程度正解できることはわかりました。
その後、50エポックの学習を行いましたが、中間層が存在しないため正解率が90%程度までしか向上しませんでした。

アンサンブル学習

活性化関数によって、正解するデータが異なるのではないかと考えました。代表的な5つの活性化関数で学習後の$ y $の平均値から正解率を求めてみます。

# 学習データ平均
y_trainm = (y_train1 + y_train2 + y_train3 + y_train4 + y_train5)/5
train_ratem = accuracy_rate(y_trainm, t_train)
# テストデータ平均
y_testm = (y_test1 + y_test2 + y_test3 + y_test4 + y_test5)/5
test_ratem = accuracy_rate(y_testm, t_test)
print("学習データ正解率 = " + str(train_ratem) + " テストデータ正解率 = " + str(test_ratem))

テスト正解率が98%を超えました。もう少し良い結果となるかと考えていましたが、やはり、判断しにくい画像はどの関数でも判断できないのでしょう。結果が少し悪かったsigmoidを除くともう少し向上しました。

活性化関数	学習正解	テスト正解
relu	100.00	97.85
sigmoid	98.35	97.38
tanh	100.00	97.89
softplus	99.87	97.68
softsign	99.97	97.79
平均	99.99	98.10
平均(sigmoidを除く)	100.00	98.18

次に、ReLU関数で、重みの初期値のシード値を変更しながら5回実行し、平均で正解率を求めます。

ReLU	学習正解	テスト正解
1	100.00	97.85
2	100.00	97.93
3	100.00	97.96
4	100.00	97.80
5	100.00	97.86
平均	100.00	98.20

判明したこと

複数の学習を組み合わせると正解率が向上する。

まとめ

今回試したパラメータのうちテスト正解率が最高のケースと正解率が97%に達したエポック数が最小のケースを表にしました。

パラメータ	正解率最高	%	97%エポック数最小	エポック
階層数2	1000	98.20	500	6
階層数3	500×100	98.30	500×50	4
多階層	6階層	98.07	6	4
重みの初期値	lecun_normal	98.07	glorot_normal	5
バイアスの初期値	0～1(乱数)	98.05	-0.1～0.1(乱数)	5
正規化関数	z_score	97.96	z_score	4
活性化関数(中間層)	tanh	97.89	relu	5
活性化関数(Leaky ReLU)	0.2	97.94	0.0	5
活性化関数(その他)	絶対値	97.82	２乗	13
活性化監修(出力層)	ソフトマックス	97.85	ソフトマックス	5
損失関数	交差エントロピー	97.85	交差エントロピー	5
学習率	0.1	97.85	0.1	5
バッチサイズ	50	97.90	10	4

テスト正解率が最高のケース(太字)と正解率が97%に達したエポック数が最小のケース(斜字)で実行しています。

パラメータ	学習正解	テスト正解	TS90%	TS95%	TS97%	TR100%	時間
正解率最高	100.00	98.16	1	1	4	13	19:18
正解率最高(relu)	100.00	98.36	1	1	2	12	18:42
97%最高	99.52	97.52	1	1	3	-	56:34
97%最高(0.01)	100.00	98.26	1	1	2	13	56:08

最高のパラメータの組み合わせでよくなるわけではないですね。正解率最高では、活性化関数をreluに変更して実行してみました。正解率97%のエポック数最小では、学習率を0.01に変更しました。結果として、学習率、エポック数とも一番よかった正解率最高(relu)の正解率グラフを描いてみます。

10エポックで、正解率98.35%に達しました。50エポックでも98.36%だったことから、10エポック実行のみで十分と言えます。

判明したこと

学習データの正解率が100%となり、テストデータの正解率が頭打ちになる。
学習率が正解率の収束には重要
ノード数を増やせば学習率が向上するが、マシンリソースを消費する。
階層を増やすと学習率の収束が早くなる。

今後の課題

学習データにフィットし、学習が進まなくなる。

先ほどのパラメータの学習データ正解率、誤差を拡大してみます。
10エポックで学習データの正解率がほぼ100%になり、誤差がほぼ0になります。学習データにフィットし、これ以上学習が進まなくなります。学習データと同時にテストデータの正解率を向上する方法を考える必要があります。

正解率の収束を早める。

次にテストデータ正解率、誤差を拡大してみます。
よく見ると正解率は、ジグザグに変化しています。重みの学習率が大きいため行ったり来たりしてるのか？また、誤差は徐々に増えています。収束の方向に向かっていません。

(参考)
プログラム全体です。前回までの関数を一部整理しています。

import gzip
import numpy as np
# MNIST読み込み
def load_mnist( mnist_path ) :
    return _load_image(mnist_path + 'train-images-idx3-ubyte.gz'), \
           _load_label(mnist_path + 'train-labels-idx1-ubyte.gz'), \
           _load_image(mnist_path + 't10k-images-idx3-ubyte.gz'), \
           _load_label(mnist_path + 't10k-labels-idx1-ubyte.gz')
def _load_image( image_path ) :
    # 画像データの読み込み
    with gzip.open(image_path, 'rb') as f:
        buffer = f.read()
    size = np.frombuffer(buffer, np.dtype('>i4'), 1, offset=4)
    rows = np.frombuffer(buffer, np.dtype('>i4'), 1, offset=8)
    columns = np.frombuffer(buffer, np.dtype('>i4'), 1, offset=12)
    data = np.frombuffer(buffer, np.uint8, offset=16)
    image = np.reshape(data, (size[0], rows[0]*columns[0]))
    image = image.astype(np.float32)
    return image
def _load_label( label_path ) :
    # 正解データ読み込み
    with gzip.open(label_path, 'rb') as f:
        buffer = f.read()
    size = np.frombuffer(buffer, np.dtype('>i4'), 1, offset=4)
    data = np.frombuffer(buffer, np.uint8, offset=8)
    label = np.zeros((size[0], 10))
    for i in range(size[0]):
        label[i, data[i]] = 1
    return label

# 正規化関数
def min_max(x, x_min=None, x_max=None, axis=None):
    if x_min is None:
        x_min = np.min(x, axis=axis, keepdims=True) # 最小値を求める
    if x_max is None:
        x_max = np.max(x, axis=axis, keepdims=True) # 最大値を求める
    return (x-x_min)/np.maximum((x_max-x_min),1e-7), x_min, x_max

def z_score(x, x_mean=None, x_std=None, axis=None):
    if x_mean is None:
        x_mean = np.mean(x, axis=axis, keepdims=True) # 平均値を求める
    if x_std is None:
        x_std  = np.std(x, axis=axis, keepdims=True)  # 標準偏差を求める
    return (x-x_mean)/np.maximum(x_std,1e-7), x_mean, x_std

# affine変換
def affine(z, W, b):
    return np.dot(z, W) + b

def affine_back(dx, z, W, b):
    dz = np.dot(dx, W.T)                # zの勾配は、今までの勾配と重みを掛けた値
    dW = np.dot(z.T, dx)                # 重みの勾配は、zに今までの勾配を掛けた値
    size = 1
    if z.ndim == 2:
        size = z.shape[0]
    db = np.dot(np.ones(size).T, dx)    # バイアスの勾配は、今までの勾配の値
    return dz, dW, db

# 活性化関数（中間）
def sigmoid(x):
    return 1/(1 + np.exp(-x))

def sigmoid_back(z):
    return z - z**2

def tanh(x):
    return np.tanh(x)

def tanh_back(x):
    return 1/np.cosh(x)**2

def relu(x):
    return np.maximum(0, x)

def relu_back(x):
    return np.where(x > 0, 1, 0)

def leaky_relu(x, alpha):
    return np.maximum(x * alpha, x)

def leaky_relu_back(x, alpha):
    return np.where(x > 0, 1, alpha)

def identity(x):
    return x

def identity_back(x):
    return 1

def softplus(x):
    return np.log(1+np.exp(x))

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

def softsign(x):
    return x/(1+np.absolute(x))

def softsign_back(x):
    return 1/(1+np.absolute(x))**2

def step(x):
    return np.where(x > 0, 1, 0)

def step_back(x):
    return x*0

# 活性化関数（出力）
def softmax(x):
    x = x.T
    max_x = np.max(x, axis=0)
    exp_x = np.exp(x - max_x)
    sum_exp_x = np.sum(exp_x, axis=0)
    y = exp_x/sum_exp_x
    return y.T

# 損失関数
def mean_squared_error(y, t):
    size = 1
    if y.ndim == 2:
        size = y.shape[0]
    return 0.5 * np.sum((y-t)**2)/size
def cross_entropy_error(y, t):
    size = 1
    if y.ndim == 2:
        size = y.shape[0]
    return -np.sum(t * np.log(np.maximum(y,1e-7)))/size
    #return -np.sum(t * np.log(y))/size

# 活性化関数(出力)＋損失関数勾配
def softmax_cross_entropy_error_back(y, x, t):
    size = 1
    if y.ndim == 2:
        size = y.shape[0]
    return (y - t)/size

# 重みの初期化
def lecun_normal(d_1, d):
    std = 1/np.sqrt(d_1)
    return np.random.normal(0, std, (d_1, d))
def lecun_uniform(d_1, d):
    min = -np.sqrt(3/d_1)
    max = np.sqrt(3/d_1)
    return np.random.uniform(min, max, (d_1, d))

def glorot_normal(d_1, d):
    std = np.sqrt(2/(d_1+d))
    return np.random.normal(0, std, (d_1, d))
def glorot_uniform(d_1, d):
    min = -np.sqrt(6/(d_1+d))
    max = np.sqrt(6/(d_1+d))
    return np.random.uniform(min, max, (d_1, d))

def he_normal(d_1, d):
    std = np.sqrt(2/d_1)
    return np.random.normal(0, std, (d_1, d))
def he_uniform(d_1, d):
    min = -np.sqrt(6/d_1)
    max = np.sqrt(6/d_1)
    return np.random.uniform(min, max, (d_1, d))

def random_w(d_1, d, min=0, max=1):
    return np.random.rand(d_1,d) *  (max-min) + min
def random_b(d, min=0, max=1):
    return np.random.rand(d) *  (max-min) + min

# 正解率
def accuracy_rate(y, t):
    max_y = np.argmax(y, axis=1)
    max_t = np.argmax(t, axis=1)
    return np.sum(max_y == max_t)/y.shape[0]

# 順伝播
def propagation(layer, x, W, b, param):
    u = {}
    z = {}
    for i in range(layer+1):
        # 入力層
        if i == 0:
            z[i] = x
        # 中間層
        elif i < layer:
            u[i] = affine(z[i-1], W[i], b[i])
            z[i] = middle_func(u[i], param)
        # 出力層
        else:
            u[i] = affine(z[i-1], W[i], b[i])
            y = output_func(u[i])
    return u, z, y

# 逆伝播
def back_propagation(layer, u, z, y, t, W, b, param):
    du = {}
    dz = {}
    dW = {}
    db = {}
    for i in range(layer, 0, -1):
        # 出力層
        if i == layer:
            du[i] = output_error_back(y, u[i], t)
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
        # 中間層
        else:
            du[i] = middle_back(dz[i], u[i], z[i], param)
            dz[i-1], dW[i], db[i] = affine_back(du[i], z[i-1], W[i], b[i])
    return du, dz, dW, db

# 活性化関数(中間)勾配
def middle_back(dz, x, z, param):
    du = middle_back_func(x, z, param)
    return dz * du

# 活性化関数(出力)＋損失関数勾配
def output_error_back(y, x, t):
    return output_error_back_func(y, x, t)

# 活性化関数(中間)定義
def middle_func(x, param):
    return relu(x)
def middle_back_func(x, z, param):
    return relu_back(x)

# 活性化関数(出力)定義
def output_func(x):
    return softmax(x)

# 損失関数定義
def error_func(y, t):
    return cross_entropy_error(y, t)

# 活性化関数(出力)＋損失関数勾配定義
def output_error_back_func(y, x, t):
    return softmax_cross_entropy_error_back(y, x, t)

import numpy as np
import time

# ノード数
d = [784, 100, 50, 10]
# 階層数
layer = len(d) - 1
# 重み、バイアスの初期化
W = {}
b = {}
for i in range(layer):
    W[i+1] = he_normal(d[i],d[i+1])
    b[i+1] = np.zeros(d[i+1])

# 学習率、バッチサイズ、エポックの設定
eta = 0.1
batch_size = 100
epoch = 50

# パラメータ設定
param = {}
#param["alpha"] = 0.2

# MNISTデータ読み込み
x_train, t_train, x_test, t_test = load_mnist('c:\\mnist\\')
# 入力データの正規化
nx_train, train_min, train_max = min_max(x_train)
nx_test,  test_min,  test_max  = min_max(x_test, train_min, train_max)

# 正解率、誤差初期化
train_rate = np.zeros(epoch+1)
test_rate = np.zeros(epoch+1)
train_err = np.zeros(epoch+1)
test_err = np.zeros(epoch+1)
# 実行(学習データ)
u_train, z_train, y_train = propagation(layer, nx_train, W, b, param)
train_rate[0] = accuracy_rate(y_train, t_train)
train_err[0] = error_func(y_train, t_train)
# 実行(テストデータ)
u_test, z_test, y_test = propagation(layer, nx_test, W, b, param)
test_rate[0] = accuracy_rate(y_test, t_test)
test_err[0] = error_func(y_test, t_test)
# 正解率、誤差表示
print(" 学習データ正解率 = " + str(train_rate[0]) + " テストデータ正解率 = " + str(test_rate[0]) +
      " 学習データ誤差 = " + str(train_err[0]) + " テストデータ誤差 = " + str(test_err[0]))

# 開始時刻設定
start_time = time.time()

for j in range(epoch):
    # データのシャッフル(正解データも同期してシャフルする必要があるため一度、結合し分離)
    nx_t = np.concatenate([nx_train, t_train], axis=1)
    np.random.shuffle(nx_t)
    nx, t = np.split(nx_t, [nx_train.shape[1]], axis=1)

    for i in range(0, nx.shape[0], batch_size):
        # 実行
        u, z, y = propagation(layer, nx[i:i+batch_size], W, b, param)
        # 勾配を計算
        du, dz, dW, db = back_propagation(layer, u, z, y, t[i:i+batch_size], W, b, param)
        
        # 重み、バイアスの調整
        for k in range(1, layer+1):
            W[k] = W[k] - eta*dW[k]
            b[k] = b[k] - eta*db[k]

    # 重み、バイアス調整後の実行(学習データ)
    u_train, z_train, y_train = propagation(layer, nx_train, W, b, param)
    train_rate[j+1] = accuracy_rate(y_train, t_train)
    train_err[j+1] = error_func(y_train, t_train)
    # 重み、バイアス調整後の実行(テストデータ)
    u_test, z_test, y_test = propagation(layer, nx_test, W, b, param)
    test_rate[j+1] = accuracy_rate(y_test, t_test)
    test_err[j+1] = error_func(y_test, t_test)
    # 正解率、誤差表示
    print(str(j+1) + " 学習データ正解率 = " + str(train_rate[j+1]) + " テストデータ正解率 = " + str(test_rate[j+1]) +
         " 学習データ誤差 = " + str(train_err[j+1]) + " テストデータ誤差 = " + str(test_err[j+1]))
    
# 終了時刻設定
end_time = time.time()
total_time = end_time - start_time
print("所要時間 = " +str(int(total_time/60))+" 分 "+str(int(total_time%60)) + " 秒")