ラビットチャレンジ レポート】深層学習 day1&day2

#0. はじめに
本記事は、Study-AIが提供するラビット・チャレンジ（JDLA認定プログラム）のレポート記事である。今回は深層学習（day1&day2）に関するレポートを記述する。

###参考書籍

1. これなら分かる深層学習入門, 瀧雅人, 講談社
2. 深層学習, 岡谷貴之, 講談社
3. ディープラーニングG検定公式テキスト, 日本ディープラーニング協会監修
4. Pytorchによる発展ディープラーニング, 小川雄太郎, マイナビ

#1. 入力層〜中間層

概要

●入力層には説明変数$\boldsymbol{x}$が入る。
●中間層（または隠れ層とも呼ばれる）の出力$f(\boldsymbol{x})$は、入力層からきた説明関数に重み$\boldsymbol{W}$を掛け合わせ、バイアス$\boldsymbol{b}$を足して活性化関数$f$により返還され、次の層へと情報が渡っていく。
$$f(\boldsymbol{x})=\boldsymbol{w}^T\boldsymbol{x}+\boldsymbol{b}$$
基本的にニューラルネットワークは、上式の繰り返しにより構築されているため、上式は重要なものとなる。

実装演習

4つのノードを持つ入力層から、3つのノードを持つ中間層への出力を実装する。式$f(\boldsymbol{x})=\boldsymbol{w}^T\boldsymbol{x}+\boldsymbol{b}$が3つ並んだイメージであるが、ベクトルであった$\boldsymbol{w}$を行列$\boldsymbol{W}$にしていることがポイント。

# 順伝播（単層・複数ユニット）

# 重み
W = np.array([
    [0.1, 0.2, 0.3], 
    [0.2, 0.3, 0.4], 
    [0.3, 0.4, 0.5],
    [0.4, 0.5, 0.6]
])

print_vec("重み", W)

# バイアス
b = np.array([0.1, 0.2, 0.3])
print_vec("バイアス", b)

# 入力値
x = np.array([1.0, 5.0, 2.0, -1.0])
print_vec("入力", x)


#  総入力
u = np.dot(x, W) + b
print_vec("総入力", u)

print_vecは、オシャレに出力してくれるように定義した独自関数である。結果は以下の通りである。

###確認テスト
#####問題1
(1).ディープラーニングは結局何をやろうとしているのか2行以内で述べよ。
明示的なプログラムの代わりに多数の中間層を持つニューラルネットワークを用いて、入力値から目的とする出力値に変換する数学モデルを構築すること。

(2)ニューラルネットワークは、どの値の最適化が最終目的か？
重みとバイアス

#####問題2
次のニューラルネットワークをかけ。
入力層：2ノード1層
中間層：3ノード2層
出力層：1ノード1層

#2. 活性化関数
###概要
●活性化関数は入力信号の総和$\boldsymbol{w}^T\boldsymbol{x}(+\boldsymbol{b})$を出力信号に変換する関数$f$のこと。
●ニューラルネットワークでは、活性化関数には非線形関数を用いないと意味がない。
　→線形関数を用いると、ニューラルネットワーク出そうを深くする意味がなくなってしまう（いつまでたっても線形だから）。
●活性化関数
・シグモイド関数
$$\sigma(x)=\dfrac{1}{1+\mathrm{e}^{-ax}}$$
・ステップ関数

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

・ReLU関数

f(x) = \left\{
\begin{array}{ll}
ｘ & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.

###実装演習
・シグモイド関数

import numpy as np
import matplotlib.pyplot as plt
from common import functions as f

x = np.array([range(-2000, 2000)])
x = x/100
y = f.sigmoid(x)

plt.scatter(x, y, marker='.')
plt.show()

・ReLU関数

import numpy as np
import matplotlib.pyplot as plt
from common import functions as f

x = np.array([range(-2000, 2000)])
x = x/100
y = f.relu(x)

plt.scatter(x, y, marker='.')
plt.show()

###確認テスト
#####問題1
配布されたソースコードより、$z=f(u)$に該当する箇所を示せ。

z1 = functions.sigmoid(u)

#3. 出力層
###誤差関数
よく使われる誤差関数には、以下の2つがある。
・二乗誤差
$$E_n(w)=\frac{1}{2}\sum_{i=1}^N(y_i-d_i)^2=\frac{1}{2}||(y-d)||^2$$
・クロスエントロピー誤差
$$E_n(w)=-\sum_{i=1}^Nd_i\log y_i$$

###活性化関数
出力層の活性化関数は、信号の大きさ（比率）をそのままに変換できる。代表的なのはソフトマックス関数であり、下記の式で表される。
$$f(i, \boldsymbol{u})=\frac{e^{u_i}}{\sum^K_{k=1}e^{u_k}}$$
このソフトマックス関数は、$i$に関して和をとると1になるという性質を持っている。この性質から、出力層の出力は0〜1の範囲に限定され、そのまま確率として利用できる（分類問題に使われる）。

###実装演習
ここではクロスエントロピーとソフトマックスを実装する。
・クロスエントロピー

def cross_entropy_error(d, y):
    if y.ndim == 1:
        d = d.reshape(1, d.size)
        y = y.reshape(1, y.size)
    # one-hot-vectorの場合、正解ラベルのインデックスに変換
    if d.size == y.size:
        d = d.argmax(axis=1)
    batch_size = y.shape[0]
    return -np.sum(np.log(y[np.arange(batch_size), d] + 1e-7)) / batch_size

・ソフトマックス

def softmax(x):
   if x.ndim == 2:#2次元だった場合
     x = x.Tx
     x = x-np.max(x, axis=0)
     y = np.exp(x) /np.sum(np.exp(x), axis=0)
      return y.T
 x = x -np.max(x) # オーバーフロー対策
      return np.exp(x) / np.sum(np.exp(x)) 

この実装から分かるように、ソフトマックス関数は通常の活性化関数とは違う点がある。通常の活性化関数はノード1つからの値を別の1つの値に変換する1対1である。しかし、ソフトマックスは出力層のすべてのノードの出力を用いて、各ノード1つずつ出力をする多対1の出力である。

###確認テスト
#####問題1
(1)2乗誤差について、なぜ引き算ではなく2乗するのかを述べよ。
単なる引き算にしてしまうと、誤差の$\pm$により誤差和の値を正確に把握することができないため。

(2)次式の$1/2$はどのような意味を持つか？
$$E_n(\boldsymbol{w})=\frac{1}{2}\sum_i(y_i-d_i)^2$$
微分した時に乗数の2が係数として出てきて、$1/2$と相殺される。微分の計算が楽になる効果がある。

4. 勾配降下法

###概要
人工知能の学習とは、出力層から計算した誤差関数の値$E(\boldsymbol{w})$を最小化することである。しかし、パラメータ$\boldsymbol{w}$のサイズが膨大で、解析的に解くのはほぼ不可能である。そこで、数値的に極小解を求める手法である勾配降下法が用いられる。
勾配降下法には、以下の3つの手法がある。

#####勾配降下法
$$w^{(t+1)}=W^{(t)}-\epsilon\nabla E$$
もっとも一般的な勾配降下法であり、別名最急降下法やニュートン法と呼ばれる。すべてのデータを用い一度に更新する手法である。

#####確率的勾配降下法（SGD）
$$w^{(t+1)}=W^{(t)}-\epsilon\nabla E_n$$
ランダムに抽出したサンプルの誤差を用いてパラメータを更新する手法である。データが冗長な場合の計算コストの軽減ができたり、望まない局所極小解に収束するリスクを軽減できる。また、学ぶべきデータが逐次増えていく場合に対応できる、オンライン学習ができる。

#####ミニバッチ勾配降下法
$$w^{(t+1)}=W^{(t)}-\epsilon\nabla E_t, E_t=\frac{1}{N_t}\sum_{n \in D_t}E_n$$
ランダムに分割したデータの集合（ミニバッチ）$D_t$に属するサンプルの平均誤差をとる。並列計算に対応できるため、確率的勾配法のメリットを損なわず、計算機の計算資源を有効利用できる。

###実装演習
ここではごく一般的な勾配降下法（最急降下法）を実装する。

class SGD:
    def __init__(self, learning_rate=0.01):
        self.learning_rate = learning_rate
        
    def update(self, params, grad):
        for key in params.keys():
            params[key] -= self.learning_rate * grad[key] 

上記プログラムから分かるように、我々が決定しなければならないパラメータは学習率$\epsilon$（上記プログラムでいうと「self.learning_rate」）である。この値を大きくすると進みは早く最小値でない極値にはまりにくいががなかなか収束せず、小さくすると極値での収束は早いが進みが遅く小さな極値も抜けれなくなってしまう。ちょうど良い値を選択するか、ベストな策としては別の手法（Adamなど）を選択するのが良い。

###確認テスト
#####問題1
(1)下記勾配降下法に該当するソースコードを配布コードから抜き出せ。
$$w^{(t+1)}=W^{(t)}-\epsilon\nabla E$$

network[key] -= learning_rate * grad[key]

(2)下記勾配降下法に該当するソースコードを配布コードから抜き出せ。
$$\nabla E=\frac{\partial E}{\partial \boldsymbol{w}}=\left[\frac{\partial E}{\partial w_1}\cdots\frac{\partial E}{\partial w_M} \right]$$

grad = backward(x, d, z1, y)

5. 誤差逆伝播法

###概要
直接数値微分計算をすると、計算量が非常に多くなってしまう。例えば推定したいパラメータの数が$n$だった場合、微分を計算しなければならない数は$n$で、さらにニューラルネットワーク全体を$n$回計算しなければならない。パラメータ数が数千万にもなるニューラルネットワークにおいて、そのような計算を毎回していては学習が進まない。
そこで、誤差逆伝播法が考え出された。算出された誤差を、出力層側から順に微分し、前の層前の層へと伝播する手法である。この手法により、最低限の計算をするだけですみ、効率よく各パラメータの更新量を解析的に計算する手法である。

###実装演習
ここでは確率的勾配降下法を誤差逆伝播法で実装する。

# サンプルとする関数
#yの値を予想するAI

def f(x):
    y = 3 * x[0] + 2 * x[1]
    return y

# 初期設定
def init_network():
    # print("##### ネットワークの初期化 #####")
    network = {}
    nodesNum = 10
    network['W1'] = np.random.randn(2, nodesNum)
    network['W2'] = np.random.randn(nodesNum)
    network['b1'] = np.random.randn(nodesNum)
    network['b2'] = np.random.randn()

    # print_vec("重み1", network['W1'])
    # print_vec("重み2", network['W2'])
    # print_vec("バイアス1", network['b1'])
    # print_vec("バイアス2", network['b2'])

    return network

# 順伝播
def forward(network, x):
    # print("##### 順伝播開始 #####")
    
    W1, W2 = network['W1'], network['W2']
    b1, b2 = network['b1'], network['b2']
    u1 = np.dot(x, W1) + b1
    z1 = functions.relu(u1)
    
    ## 試してみよう
    #z1 = functions.sigmoid(u1)
    
    u2 = np.dot(z1, W2) + b2
    y = u2

    # print_vec("総入力1", u1)
    # print_vec("中間層出力1", z1)
    # print_vec("総入力2", u2)
    # print_vec("出力1", y)
    # print("出力合計: " + str(np.sum(y)))    
    
    return z1, y

# 誤差逆伝播
def backward(x, d, z1, y):
    # print("\n##### 誤差逆伝播開始 #####")    

    grad = {}
    
    W1, W2 = network['W1'], network['W2']
    b1, b2 = network['b1'], network['b2']

    # 出力層でのデルタ
    delta2 = functions.d_mean_squared_error(d, y)
    # b2の勾配
    grad['b2'] = np.sum(delta2, axis=0)
    # W2の勾配
    grad['W2'] = np.dot(z1.T, delta2)
    # 中間層でのデルタ
    #delta1 = np.dot(delta2, W2.T) * functions.d_relu(z1)

    ## 試してみよう
    delta1 = np.dot(delta2, W2.T) * functions.d_sigmoid(z1)

    delta1 = delta1[np.newaxis, :]
    # b1の勾配
    grad['b1'] = np.sum(delta1, axis=0)
    x = x[np.newaxis, :]
    # W1の勾配
    grad['W1'] = np.dot(x.T, delta1)
    
    # print_vec("偏微分_重み1", grad["W1"])
    # print_vec("偏微分_重み2", grad["W2"])
    # print_vec("偏微分_バイアス1", grad["b1"])
    # print_vec("偏微分_バイアス2", grad["b2"])

    return grad

# サンプルデータを作成
data_sets_size = 100000
data_sets = [0 for i in range(data_sets_size)]

for i in range(data_sets_size):
    data_sets[i] = {}
    # ランダムな値を設定
    data_sets[i]['x'] = np.random.rand(2)
    
    ## 試してみよう_入力値の設定
    # data_sets[i]['x'] = np.random.rand(2) * 10 -5 # -5〜5のランダム数値
    
    # 目標出力を設定
    data_sets[i]['d'] = f(data_sets[i]['x'])
    
losses = []
# 学習率
learning_rate = 0.07

# 抽出数
epoch = 1000

# パラメータの初期化
network = init_network()
# データのランダム抽出
random_datasets = np.random.choice(data_sets, epoch)

# 勾配降下の繰り返し
for dataset in random_datasets:
    x, d = dataset['x'], dataset['d']
    z1, y = forward(network, x)
    grad = backward(x, d, z1, y)
    # パラメータに勾配適用
    for key in ('W1', 'W2', 'b1', 'b2'):
        network[key]  -= learning_rate * grad[key]

    # 誤差
    loss = functions.mean_squared_error(d, y)
    losses.append(loss)

print("##### 結果表示 #####")    
lists = range(epoch)


plt.plot(lists, losses, '.')
# グラフの表示
plt.show()

結果は以下の通りになった。横軸はエポックで、縦軸は損失を表している。

###確認テスト
#####問題1
（1)下記式を表すソースコードを探せ。
$$\frac{\partial E}{\partial \boldsymbol{y}}$$

delta2 = functions.d_mean_squared_error(d, y)

$$\frac{\partial E}{\partial \boldsymbol{y}}\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{u}}$$

delta2 = functions.d_mean_squared_error(d, y)

$$\frac{\partial E}{\partial \boldsymbol{y}}\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{u}}\frac{\partial \boldsymbol{u}}{\partial w_{ji}^{(2)}}$$

grad['W2'] = np.dot(z1.T, delta2)

6. 勾配消失問題

###概要
勾配消失問題とは、デルタの逆伝播にまつわる深刻な問題である。順伝播においては、ノードへの入力は活性化関数により変換されたのち、次層への入力となる。活性化関数がシグモイドなら、入力と出力の関係は
$$\boldsymbol{z}^{(l)}=\sigma(\boldsymbol{W}^{(l)}\sigma(\boldsymbol{W}^{(l-1)}\cdots\sigma(\boldsymbol{W}^{(1)}\boldsymbol{z}^{(0)})\cdots))$$
という非線形な合成関数の関係にある。活性化関数がシグモイドの場合を考えると、各ユニットの出力が0から1の間の値に制限されているために、信号の大きさが伝播中に爆発する心配はない。また、ユニットへの層入力がかなり小さな値になってしまっても、シグモイドによってある程度大きな出力地へ引き延ばされるため、信号が消えてしまう心配も一般にはない。
それと比べて逆伝播は、デルタを活性化関数で変換することはなく線形な変換である。というのも、デルタの値は
$$\delta_j^{(l)}=\sum_{q, p, \cdots, k}\delta_q^{(L)}w_{qp}^{(L)}f'(u_{p}^{(L-1)})\cdots w_{lk}^{(l+2)}f'(u_{k}^{(l+1)})w_{kj}^{(l+1)}f'(u_{j}^{(l)})$$
というように1区間ごとに活性化関数の微分値$f'(u)$倍された後、そのまま下層へと伝播されるからである。したがって、もし出力層からある中間層まで$m$区間だけ逆伝播したとすると、デルタの値は出力層のものと比べておよそ$(f'(u))^m$倍されている。この指数的な因子のせいで、逆伝播中にデルタの値は急激に消失してしまう。
シグモイドの微分は1/4以下であるという性質を持っているため、$(f'(u))^m=(1/4)^m$倍以下になってしまう。例えば$m=10$であると約0.00001倍だけ減衰されてしまう。デルタの値は指数的にどんどん小さくなってしまう。そのため、入力層側へ向かうにつれて誤差関数の勾配の値もどんどん小さくなり、勾配を用いたパラメータの更新が一向に進まなくなる、つまり勾配降下法を用いても。入力層よりの領域では重みパラメータは全く更新されず、ニューラルネットの学習がうまくいかなくなってしまう。

###実装演習
下記で意図的に勾配消失を起こす。初期値はガウス関数をもとに生成し、活性化関数にはシグモイドを使っている。

import numpy as np
from common import layers
from collections import OrderedDict
from common import functions
from data.mnist import load_mnist
import matplotlib.pyplot as plt

# mnistをロード
(x_train, d_train), (x_test, d_test) = load_mnist(normalize=True, one_hot_label=True)
train_size = len(x_train)

print("データ読み込み完了")

# 重み初期値補正係数
wieght_init = 0.01
#入力層サイズ
input_layer_size = 784
#中間層サイズ
hidden_layer_1_size = 40
hidden_layer_2_size = 20

#出力層サイズ
output_layer_size = 10
# 繰り返し数
iters_num = 2000
# ミニバッチサイズ
batch_size = 100
# 学習率
learning_rate = 0.1
# 描写頻度
plot_interval=10

# 初期設定
def init_network():
    network = {} 
    network['W1'] = wieght_init * np.random.randn(input_layer_size, hidden_layer_1_size)
    network['W2'] = wieght_init * np.random.randn(hidden_layer_1_size, hidden_layer_2_size)
    network['W3'] = wieght_init * np.random.randn(hidden_layer_2_size, output_layer_size)

    network['b1'] = np.zeros(hidden_layer_1_size)
    network['b2'] = np.zeros(hidden_layer_2_size)
    network['b3'] = np.zeros(output_layer_size)

    return network

# 順伝播
def forward(network, x):
    W1, W2, W3 = network['W1'], network['W2'], network['W3']
    b1, b2, b3 = network['b1'], network['b2'], network['b3']
    hidden_f = functions.sigmoid
    
    u1 =  np.dot(x, W1) + b1
    z1 = hidden_f(u1)
    u2 =  np.dot(z1, W2) + b2
    z2 = hidden_f(u2)
    u3 =  np.dot(z2, W3) + b3
    y = functions.softmax(u3)
 
    return z1, z2, y

# 誤差逆伝播
def backward(x, d, z1, z2, y):
    grad = {}
    
    W1, W2, W3 = network['W1'], network['W2'], network['W3']
    b1, b2, b3 = network['b1'], network['b2'], network['b3']
    hidden_d_f = functions.d_sigmoid
    last_d_f = functions.d_softmax_with_loss
    
    
    # 出力層でのデルタ
    delta3 = last_d_f(d, y)
    # b3の勾配
    grad['b3'] = np.sum(delta3, axis=0)
    # W3の勾配
    grad['W3'] = np.dot(z2.T, delta3)
    # 2層でのデルタ
    delta2 = np.dot(delta3, W3.T) * hidden_d_f(z2)
    # b2の勾配
    grad['b2'] = np.sum(delta2, axis=0)
    # W2の勾配
    grad['W2'] = np.dot(z1.T, delta2)
    # 1層でのデルタ
    delta1 = np.dot(delta2, W2.T) * hidden_d_f(z1)
    # b1の勾配
    grad['b1'] = np.sum(delta1, axis=0)
    # W1の勾配
    grad['W1'] = np.dot(x.T, delta1)

    return grad

# パラメータの初期化
network = init_network()

accuracies_train = []
accuracies_test = []

# 正答率
def accuracy(x, d):
    z1, z2, y = forward(network, x)
    y = np.argmax(y, axis=1)
    if d.ndim != 1 : d = np.argmax(d, axis=1)
    accuracy = np.sum(y == d) / float(x.shape[0])
    return accuracy

for i in range(iters_num):
    # ランダムにバッチを取得    
    batch_mask = np.random.choice(train_size, batch_size)
    # ミニバッチに対応する教師訓練画像データを取得    
    x_batch = x_train[batch_mask]
    # ミニバッチに対応する訓練正解ラベルデータを取得する
    d_batch = d_train[batch_mask]


    
    z1, z2, y = forward(network, x_batch)
    grad = backward(x_batch, d_batch, z1, z2, y)

    if (i+1)%plot_interval==0:
        accr_test = accuracy(x_test, d_test)
        accuracies_test.append(accr_test)
        
        accr_train = accuracy(x_batch, d_batch)
        accuracies_train.append(accr_train)

        print('Generation: ' + str(i+1) + '. 正答率(トレーニング) = ' + str(accr_train))
        print('                : ' + str(i+1) + '. 正答率(テスト) = ' + str(accr_test))

    # パラメータに勾配適用
    for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
        network[key]  -= learning_rate * grad[key]


lists = range(0, iters_num, plot_interval)
plt.plot(lists, accuracies_train, label="training set")
plt.plot(lists, accuracies_test,  label="test set")
plt.legend(loc="lower right")
plt.title("accuracy")
plt.xlabel("count")
plt.ylabel("accuracy")
plt.ylim(0, 1.0)
# グラフの表示
plt.show()

結果は以下の通りである。

見て分かるように、どんなに学習しても一向に正解率が上がらない。
では、初期値の生成方法はそのままガウス関数で、活性化関数をReLUにしてみたらどうか？結果はかきになる（実装は省略）

勾配消失は起きず、学習するにつれてしっかりと正解率が上がっていることがわかる。
このように、勾配消失と活性化関数は密接に関わっていることがわかる。

###確認テスト
#####問題1
(1)連鎖律の原理を使い、$z=t^2$、$t=x+y$の時の$dz/dx$を求めよ。

\frac{dz}{dx}=\frac{dz}{dt}\frac{dt}{dz}=2t\times1=2(x+y)

#####問題2
(1)シグモイド関数を微分した時、入力値が0の時に最大値をとる、その値を答えよ。
0.25

7. 学習率最適化手法

###概要
#####モメンタム（Momentum）
モーメンタムとは運動量という意味の言葉で、物理に関係がある。これは勾配降下法の収束性能を向上させる方法の一つであり、重みの修正量に前回の重みの修正量のいくらかを加算する方法で、ミニバッチ$t-1$に対する重みの修正量を$\Delta\boldsymbol{w}^{(t-1)}\equiv\boldsymbol{w}^{(t)}-\boldsymbol{w}^{(t-1)}$と書くと、ミニバッチ$t$に対する更新を

\begin{align}
\boldsymbol{w}^{(t+1)}&=\boldsymbol{w}^{(t)}-\epsilon\nabla E_t+\mu\Delta \boldsymbol{w}^{(t-1)}\\
\Delta\boldsymbol{w}^{(t+1)}&=\mu\Delta\boldsymbol{w}^{(t-1)}-\epsilon\nabla E_t
\end{align}

と定めるのである。$\mu$は加算の割合を制御するハイパーパラメータであり、通常$\mu=0.5\sim0.9$程度の範囲から選ぶ。
誤差関数が深い谷状の形状を持ち、かつその谷底にあまり高低差がないときには、勾配降下法は非常に効率が悪くなることが知られており。谷が深いので谷底を少しでも外れた点では谷と直交する方向に大きな勾配が生じる。その結果重みは毎度、谷と直交する方向に修正され、その経路はジグザグになり谷底を効率よく探索することができない。しかしこの手法では、修正量が過去の修正量の重み付き平均になるため、谷と直交する方向の修正量は平均化されてなくなり、重み修正の経路はジグザグではなく谷底を谷の方向に沿って一直線に辿ることになる。

#####学習係数を変化させる勾配降下法（AdaGrad）
勾配降下法では、パラメータの更新量の大きさは学習係数によって変わる。この学習係数をどう決めるかで学習の成否が左右されるため、極めて重要なパラメータである。学習係数の決め方を論じた研究はいくつもあり、また自動的に決定する方法も提案されているが、今でも手動で（試行錯誤的に）値を選ぶことが一般的であると言える。
そのような学習係数を決める際によく用いられる方法が、学習が進むにつれて学習係数を小さくするという方法である。最初は大きく学習し、次第に小さく学習するという手法でニューラルネットワークの学習では実際によく使われる。
学習係数を徐々に下げていくというアイデアは、パラメータ「全体」の学習係数の値を一括して下げることに相当している。これをさらに発展させたのがAdaGradである。AdaGradは「一つひとつ」のパラメータに対して、オーダーメイドの値を設定することが特徴である。この手法では、誤差関数の勾配を$g_{t}\equiv\nabla E_t$と書き、この成分を$g_{t,i}$と書くと、普通の勾配降下法では更新量の$i$成分は$-\epsilon g_{t,i}$だが、AdaGradではこれを

\begin{align}
-\dfrac{\epsilon}{\sqrt{\sum_{t'=1}^t g_{t',i}^2}}g_{t,i}
\end{align}

とする。直感的には、大きく変化する成分の更新を早くに制御していく、ということになる。AdaGradの欠点としては、学習の初期に勾配が大きいとすぐさま更新量が小さくなってしまうことである。そのため、適切な程度にまで$\epsilon$を大きく選ぶ必要があり、学習係数の選び方に敏感で使いにくい手法となっている。

#####AdaGradの改善法：RMSprop
AdaGradでは、過去のすべての勾配の情報を集積してしまっていたため、非常に扱いにくい面があった。しかし、RMSpropという方法では、十分過去の勾配の情報を指数関数的な元帥因子によって消滅させられるように、二乗和ではなく指数的な移動平均$v_{t,i}$から決まるroot mean square（RMS）を用いることにした。

\begin{align}
v_{t,i}=&\rho\, v_{t-1,i}+(1-\rho)(\nabla E(\boldsymbol{w}^{(t)})_i)^2\\
\boldsymbol{w}^{(t+1)}=&\boldsymbol{w}^{(t)}-\dfrac{\epsilon}{\sqrt{v_{t,i}+\eta}}\nabla E(\boldsymbol{w}^{(t)})_i
\end{align}

RMSpropでは最近の勾配の履歴のみが影響するため、更新量が完全に消えることがないという強みを持つ。

#####RMSpropの改良法：Adam
最後に、RMSpropの改良法を説明する。Adamでは、分母にある勾配のRMS（$\sqrt{v_{t,i}+\eta}$）のみならず、勾配自身も指数的な移動平均による推定値で置き換える手法である。これはRMSpropの勾配部分に、指数的減衰を含むモメンタムを適用したようなものだが、実際はAdamはかなり手が込んでいる。
まず勾配とその2乗の指数的な移動平均を定義する。

\begin{align}
m_{t,i}=&\rho_1\, m_{t-1,i}+(1-\rho_1)\nabla E(\boldsymbol{w}^{(t)})_i\\
v_{t,i}=&\rho_2\, v_{t-1,i}+(1-\rho_2)(\nabla E(\boldsymbol{w}^{(t)})_i)^2
\end{align}

ただし初期値は$m_{0,i}=v_{0,i}=0$である。これは一見すると勾配の1次モーメントと2次モーメントの良い推定量のように見えるが、実はバイアスを持っている。というのも初期値は0にとってしまうので、更新の初期はモーメントの推定量が0の方に偏ってしまう。
そこでバイアスを補正し、できるだけ不偏推定量に近づくようにしたのが以下の式である（めんどくさいので詳細は省く）。するとバイアス補正したモーメントの推定値は

\begin{align}
\hat{m}_{t,i}=\dfrac{m_{t,i}}{(1-\rho_1^t)},\hspace{3mm}\hat{v}_{t,i}=\dfrac{v_{t,i}}{(1-\rho_2^t)}
\end{align}

となる。Adamはこれらの推定値を用いた勾配降下法である。

\begin{align}
\boldsymbol{w}^{(t+1)}=&\boldsymbol{w}^{(t)}-\epsilon\dfrac{\hat{m}_{t,i}}{\sqrt{\hat{v}_{t,i}+\eta}}
\end{align}

Adamの原論文におけるパラメータの推奨値は
$$\epsilon=0.001,\hspace{2mm}\rho_1=0.9,\hspace{2mm}\rho_2=0.999,\hspace{2mm}\eta=10^{-8}$$
である。さまざまな深層学習のフレームワークでも、基本的にはこの推奨値が用いられている。

###実装演習
上記の数式を実装していく。
#####モメンタム

class Momentum:
    def __init__(self, learning_rate=0.01, momentum=0.9):
        self.learning_rate = learning_rate
        self.momentum = momentum
        self.v = None
        
    def update(self, params, grad):
        if self.v is None:
            self.v = {}
            for key, val in params.items():                                
                self.v[key] = np.zeros_like(val)
                
        for key in params.keys():
            self.v[key] = self.momentum * self.v[key] - self.learning_rate * grad[key] 
            params[key] += self.v[key]

#####AdaGrad

class AdaGrad:
    def __init__(self, learning_rate=0.01):
        self.learning_rate = learning_rate
        self.h = None
        
    def update(self, params, grad):
        if self.h is None:
            self.h = {}
            for key, val in params.items():
                self.h[key] = np.zeros_like(val)
            
        for key in params.keys():
            self.h[key] += grad[key] * grad[key]
            params[key] -= self.learning_rate * grad[key] / (np.sqrt(self.h[key]) + 1e-7)

#####RMSprop

class RMSprop:
    def __init__(self, learning_rate=0.01, decay_rate = 0.99):
        self.learning_rate = learning_rate
        self.decay_rate = decay_rate
        self.h = None
        
    def update(self, params, grad):
        if self.h is None:
            self.h = {}
            for key, val in params.items():
                self.h[key] = np.zeros_like(val)
            
        for key in params.keys():
            self.h[key] *= self.decay_rate
            self.h[key] += (1 - self.decay_rate) * grad[key] * grad[key]
            params[key] -= self.learning_rate * grad[key] / (np.sqrt(self.h[key]) + 1e-7)

#####Adam

class Adam:
    def __init__(self, learning_rate=0.001, beta1=0.9, beta2=0.999):
        self.learning_rate = learning_rate
        self.beta1 = beta1
        self.beta2 = beta2
        self.iter = 0
        self.m = None
        self.v = None
        
    def update(self, params, grad):
        if self.m is None:
            self.m, self.v = {}, {}
            for key, val in params.items():
                self.m[key] = np.zeros_like(val)
                self.v[key] = np.zeros_like(val)
        
        self.iter += 1
        lr_t  = self.learning_rate * np.sqrt(1.0 - self.beta2 ** self.iter) / (1.0 - self.beta1 ** self.iter)         
        
        for key in params.keys():
            self.m[key] += (1 - self.beta1) * (grad[key] - self.m[key])
            self.v[key] += (1 - self.beta2) * (grad[key] ** 2 - self.v[key])
            
            params[key] -= lr_t * self.m[key] / (np.sqrt(self.v[key]) + 1e-7)

###確認テスト
#####問題1
(1)モメンタム・AdaGrad・RMSPropの特徴をそれぞれ簡潔に説明せよ。
●モメンタムのメリット
・局所的最適解にはならず、大域的最適解になる
・谷間についてから最も低い位置（最適値）に行くまでの時間が早い。

●AdaGradのメリット
・勾配の緩やかな斜面に対して、最適値に近づける。

●RMSPropのメリット
・局所的最適解にはならず、大域的最適解となる。
・ハイパーパラメータの調整が必要な場合が少ない。

8.過学習

###概要
過学習とはテスト誤差と訓練誤差とで学習曲線が乖離することである。訓練誤差は下がっていくのに対し、テスト誤差は下がらずに上がってしまう現象である。原因としては以下が考えられる。
・パラメータの数が多い
・パラメータの値が適切ではない
・ノードが多い
などなど。つまりこれは、ネットワークの自由度が高いということになる。では、それを防ぐ方法としては何があるか？

#####pノルム
下記の式で表されるノルムを、pノルムという。
$$||x||_p=\left(|x_1|^p+\cdots+|x_n|^p\right)^{1/p}$$
$p=1$の時、L1正則化と呼ぶ。$p=2$の時は、L2正則化という。

#####ドロップアウト
過学習の課題は、ネットワークの自由度が高いことにより起こっていた。そこでドロップアウトを使うことで、学習中にランダムにノードを削除して適度な自由度下で学習させることができる。またメリットとして、データ量を変化させずに、異なるモデルを学習させていると解釈できる。

###実装演習
単純な過学習、L2正則化、L1正則化、dropoutを実装した。

これを見ると、それぞれに特徴があることが見て取れる。過学習は訓練の正解率が1となっているのに対してテストの正解率はそこまで上がってない。L2正則化は訓練とテストの正解率の乖離が小さくなっている。L1正則化は正解率が下がったり上がったりして、収束しづらいことがわかる。そしてdropoutは緩やかな正解率の上昇ではあるが、安定して正解率が上がっている。

###確認テスト
#####問題1
(1)下図について、L1正則化を表しているグラフはどちらか？

左がL1正則化を表している。L1正則化は絶対値の形になるので、図形にすると角がある（四角の）形になる。L2正則化は2乗になっているので、図にすると円で表されることがわかる。

9. 畳み込みニューラルネットワークの概念

###概要
#####フィルターによる計算
畳み込みニューラルネットワーク(CNN)とは、通常のニューラルネットワークとは違い、画像を扱うのに特化したニューラルネットワークである。通常のニューラルネットワークでは画像の上下左右斜めのピクセルの関係を維持して学習することはできないが（入力が一直線に並んだニューロンであるため）、CNNであればした図で表すような（重み）フィルターを用いることでその課題を解決し、画像の扱いに特化した。

下記のように、フィルターと入力画像を重ね合わせ、重なった部分同士の積をとり、それらを全て足し合わせて一個の出力にする。上記の例では以下のような式と出力になる。
$$3\times3+4\times1+4\times2+0\times8+8\times7+9\times5+0\times5+4\times4+3\times1=141$$

この操作をフィルターを1つずつずらしながら全ての入力画像の点を通るように行う。つまり下記のような操作である（計算過程は略）。

#####パディング
パディングとは、入力画像の周りに適当な値のピクセルを増やすことにより、入力画像を大きくする操作である。通常CNNでは入力画像より出力の方が小さくなるので、繰り返すと画像がかなり小さくなってしまう。その対策として、入力画像を大きくすることで出力もある程度の大きさを保てるのである。
ちなみに、埋める値はなんでも良いが、全てをゼロで埋めることをゼロパディングという。

#####ストライド
通常のCNNはフィルターを入力画像に沿って1ずつスライドさせるが、ストライドの値を変えることにより2ずつだったり3ずつだったり、フィルターの移動距離を伸ばすことができる。
しかしその場合、もちろん出力の値はどんどん小さくなっていくので注意である。

入力サイズを$W\times H$、フィルタサイズを$Fw\times Fh$、パディングを$p$、ストライドを$s$とし、畳み込み層の出力サイズを$OW\times OH$とすると、$OW$および$OH$は次の式により求められる。
$$OW=\frac{W+2p−Fw}{s}+1$$
$$OH=\frac{H+2p−Fh}{s}+1$$

#####チャンネル
今まで見てきた例では入力画像は$W\times H$の2次元画像だったが、3次元画像にしても良い。その場合はフィルターも3次元にする必要がある。
チャンネルを増やす例は、例えばRGBのカラー画像を学習するときである。R、G、Bのそれぞれを各チャンネルとして表し、CNNへの入力として扱うのである。

###実装演習
上記概要で学んだ畳み込み操作のフォーワードとバックプロパゲーションの実装は以下のようになる。

class Convolution:
    # W: フィルター, b: バイアス
    def __init__(self, W, b, stride=1, pad=0):
        self.W = W
        self.b = b
        self.stride = stride
        self.pad = pad
        
        # 中間データ（backward時に使用）
        self.x = None   
        self.col = None
        self.col_W = None
        
        # フィルター・バイアスパラメータの勾配
        self.dW = None
        self.db = None

    def forward(self, x):
        # FN: filter_number, C: channel, FH: filter_height, FW: filter_width
        FN, C, FH, FW = self.W.shape
        N, C, H, W = x.shape
        # 出力値のheight, width
        out_h = 1 + int((H + 2 * self.pad - FH) / self.stride)
        out_w = 1 + int((W + 2 * self.pad - FW) / self.stride)
        
        # xを行列に変換
        col = im2col(x, FH, FW, self.stride, self.pad)
        # フィルターをxに合わせた行列に変換
        col_W = self.W.reshape(FN, -1).T

        out = np.dot(col, col_W) + self.b
        # 計算のために変えた形式を戻す
        out = out.reshape(N, out_h, out_w, -1).transpose(0, 3, 1, 2)

        self.x = x
        self.col = col
        self.col_W = col_W

        return out

    def backward(self, dout):
        FN, C, FH, FW = self.W.shape
        dout = dout.transpose(0, 2, 3, 1).reshape(-1, FN)

        self.db = np.sum(dout, axis=0)
        self.dW = np.dot(self.col.T, dout)
        self.dW = self.dW.transpose(1, 0).reshape(FN, C, FH, FW)

        dcol = np.dot(dout, self.col_W.T)
        # dcolを画像データに変換
        dx = col2im(dcol, self.x.shape, FH, FW, self.stride, self.pad)

        return dx

上のクラスを用いて実際に画像をCNNに学習させてみると、以下の結果となった。訓練の正解率とテストの正解率の差が小さく、どちらも90%を超えていることから、かなり精度よく学習できていることがわかる。

###確認テスト
#####問題1
(1)サイズ$6\times6$の入力画像を、サイズ$2\times2$のフィルタで畳み込んだ時の出力画像のサイズを答えよ。なお、ストライドとパディングは1とする。
求める出力幅を$OW$、出力高を $OH$とすると、サイズを求める公式により以下のようになる。
$$OW=\frac{W+2p−Fw}{s}+1=\frac{6+2−2}{1}+1=7$$
$$OH=\frac{H+2p−Fh}{s}+1=\frac{6+2−2}{1}+1=7$$

10. 最新のCNN

###概要
ここでは最新のCNNとしてAlexNetを紹介する。

論文の筆頭著者Alex Krizhevskyの名前から、AlexNetと名づけられている。AlexNetは5層の畳み込み層およびプーリング層など、それに続く3層の全結合層から構成されている。現在よく使われているCNNと比べると浅いニューラルネットワークになっているが、Yann LeCunらによって1998年に初めて考案されたCNNであるLeNetと比較すると、かなり深い構造になっている。
そのため過学習も比較的起きやすくなっていた。AlexNetでは、過学習を防ぐために、サイズ4096の全結合層の出力にドロップアウトを使用している。

しかしAlexNetは、2021年現在AlexNetは全く最新と言えない状況である。なぜならAlexNetは、2012年のILSVRCにおいて圧倒的な精度を誇ったモデルだからである。AlexNet行こう、畳み込みとプーリングの繰り返しをさらに増やした、より深いネットワークのモデルが続々と登場した。AlexNetよりもさらに深いVGGやGoogLeNetというモデルは、ILSVRCの記録をさらに大きく塗り替えていった。

###実装演習
AlexNetをPytorchにて実装したコードが以下のようになる。

import torch
import torchvision
import torch.nn as nn
import torch.nn.init as init
import torch.optim as optim
import torch.nn.functional as F
import torchvision.transforms as transforms
import numpy as np
from matplotlib import pyplot as plt
 
# load CIFA-10 data
train_dataset = torchvision.datasets.CIFAR10(
      root='./data/',
      train=True, 
      transform=transforms.ToTensor(),
      download=True)
 
test_dataset = torchvision.datasets.CIFAR10(
      root='./data/',
      train=False, 
      transform=transforms.ToTensor(),
      download=True)
 
print ('train_dataset = ', len(train_dataset))
print ('test_dataset = ', len(test_dataset))
image, label = train_dataset[0]
print (image.size())
 
# set data loadser
train_loader = torch.utils.data.DataLoader(
      dataset=train_dataset,
      batch_size=64, 
      shuffle=True,
      num_workers=2)
 
test_loader = torch.utils.data.DataLoader(
      dataset=test_dataset,
      batch_size=64, 
      shuffle=False,
      num_workers=2)
 
# Alexnet
class AlexNet(nn.Module):
 
    def __init__(self, num_classes):
        super(AlexNet, self).__init__()
        self.features = nn.Sequential(
            nn.Conv2d(3, 64, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=2, stride=2),
            nn.Conv2d(64, 192, kernel_size=5, padding=2),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=2, stride=2),
            nn.Conv2d(192, 384, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.Conv2d(384, 256, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.Conv2d(256, 256, kernel_size=3, padding=1),
            nn.ReLU(inplace=True),
            nn.MaxPool2d(kernel_size=2, stride=2),
        )
        self.classifier = nn.Sequential(
            nn.Dropout(),
            nn.Linear(256 * 4 * 4, 4096),
            nn.ReLU(inplace=True),
            nn.Dropout(),
            nn.Linear(4096, 4096),
            nn.ReLU(inplace=True),
            nn.Linear(4096, num_classes),
        )
 
    def forward(self, x):
        x = self.features(x)
        x = x.view(x.size(0), 256 * 4 * 4)
        x = self.classifier(x)
        return x
 
# select device
num_classes = 10
device = 'cuda' if torch.cuda.is_available() else 'cpu'
net = AlexNet(num_classes).to(device)
 
# optimizing
criterion = nn.CrossEntropyLoss()
optimizer = optim.SGD(net.parameters(), lr=0.01, momentum=0.9, weight_decay=5e-4)
 
# training
num_epochs = 20
train_loss_list, train_acc_list, val_loss_list, val_acc_list = [], [], [], []
 
### training
for epoch in range(num_epochs):
    train_loss, train_acc, val_loss, val_acc = 0, 0, 0, 0
    
    # ====== train_mode ======
    net.train()
    for i, (images, labels) in enumerate(train_loader):
      images, labels = images.to(device), labels.to(device)      
      optimizer.zero_grad()
      outputs = net(images)
      loss = criterion(outputs, labels)
      train_loss += loss.item()
      train_acc += (outputs.max(1)[1] == labels).sum().item()
      loss.backward()
      optimizer.step()
    
    avg_train_loss = train_loss / len(train_loader.dataset)
    avg_train_acc = train_acc / len(train_loader.dataset)
    
    # ====== val_mode ======
    net.eval()
    with torch.no_grad():
      for images, labels in test_loader:
        images = images.to(device)
        labels = labels.to(device)
        outputs = net(images)
        loss = criterion(outputs, labels)
        val_loss += loss.item()
        val_acc += (outputs.max(1)[1] == labels).sum().item()
    avg_val_loss = val_loss / len(test_loader.dataset)
    avg_val_acc = val_acc / len(test_loader.dataset)
    
    print ('Epoch [{}/{}], Loss: {loss:.4f}, val_loss: {val_loss:.4f}, val_acc: {val_acc:.4f}' 
                   .format(epoch+1, num_epochs, i+1, loss=avg_train_loss, val_loss=avg_val_loss, val_acc=avg_val_acc))
    train_loss_list.append(avg_train_loss)
    train_acc_list.append(avg_train_acc)
    val_loss_list.append(avg_val_loss)
    val_acc_list.append(avg_val_acc)
    
 
# plot graph
plt.figure()
plt.plot(range(num_epochs), train_loss_list, color='blue', linestyle='-', label='train_loss')
plt.plot(range(num_epochs), val_loss_list, color='green', linestyle='--', label='val_loss')
plt.legend()
plt.xlabel('epoch')
plt.ylabel('loss')
plt.title('Training and validation loss')
plt.grid()
 
plt.figure()
plt.plot(range(num_epochs), train_acc_list, color='blue', linestyle='-', label='train_acc')
plt.plot(range(num_epochs), val_acc_list, color='green', linestyle='--', label='val_acc')
plt.legend()
plt.xlabel('epoch')
plt.ylabel('acc')
plt.title('Training and validation accuracy')
plt.grid()