【ラビットチャレンジ（E資格）】深層学習（day2）

#はじめに
2021/2/19・20に実施される日本ディープラーニング協会(JDLA)Ｅ資格合格を目指して、ラビットチャレンジを受講した際の学習記録です。

###科目一覧
応用数学
機械学習
深層学習（day1）
深層学習（day2）
深層学習（day3）
深層学習（day4）

#Section1：勾配消失問題
誤差逆伝播法が下位層に進んでいくにつれて、勾配がどんどん緩やかになっていく。
そのため、勾配降下法による更新では下位層のパラメータはほとんど変わらず、訓練は最適値に収束しなくなる。

###勾配消失の解決方法

• 活性化関数の選択

  • ReLU関数
    $$ f(x) = \left\{ \begin{array} \
    x & (x > 0) \
    0 & (x \leq 0) \
    \end{array} \right. $$
    勾配消失問題の回避とスパース化に貢献することで良い成果をもたらしている。

• 重みの初期値設定

  • Xavier：前の層のノード数がnであるとき、 重みの要素に　$\sqrt{\frac{1}{n}}$ を乗算した値。活性化関数はReLu、シグモイド（ロジスティック）関数、双曲線正接関数（tanh）。
  • He：前の層のノード数がnであるとき、 重みの要素に　$\sqrt{\frac{2}{n}}$ を乗算した値。活性化関数はReLu。

重みの初期値に0を設定すると、どのような問題が発生するか？
→　全ての値が同じ値で伝わるため。パラメータのチューニングが行われなくなる。

• バッチ正規化
  ミニバッチ単位で入力値のデータの偏りを抑制する手法。
  中間層出力を正規化する処理を孕んだ層を加えることで、出力が常に平均0、分散1の分布に従うように強制する。
  計算の高速化、勾配消失が起きづらくなるというメリットがある。

ミニバッチの平均と分散は
$$\mu_t=\frac{1}{N_t}\sum_{i=1}^{N_t}x_{ni}, \quad \sigma_t^2=\frac{1}{N_t}\sum_{i=1}^{N_t}(x_{ni}-\mu_t)^2$$
と表され、出力を正規化すると
$$\hat x_{ni}=\frac{x_{ni}-\mu_t}{\sqrt{\sigma_t^2-\theta}}$$
となる。
この正規化された出力を学習可能なスケーリングパラメータ $\gamma$ 、シフトパラメータ $\beta$ で線形変換する。
$$y_{ni}=\gamma x_{ni}+\beta$$

活性化関数、重みの初期値、バッチ正規化の有無を変更した場合の勾配消失への寄与

import sys, os
sys.path.append(os.pardir)  # 親ディレクトリのファイルをインポートするための設定
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
from common import optimizer

class MultiLayerNet:
    '''
    input_size: 入力層のノード数
    hidden_size_list: 隠れ層のノード数のリスト
    output_size: 出力層のノード数
    activation: 活性化関数
    weight_init_std: 重みの初期化方法
    weight_decay_lambda: L2正則化の強さ
    use_dropout: ドロップアウトの有無
    dropout_ratio: ドロップアウト率
    use_batchnorm: バッチ正規化の有無
    '''
    def __init__(self, input_size, hidden_size_list, output_size, activation='relu', weight_init_std='relu', weight_decay_lambda=0,
                 use_dropout = False, dropout_ratio = 0.5, use_batchnorm=False):
        self.input_size = input_size
        self.output_size = output_size
        self.hidden_size_list = hidden_size_list
        self.hidden_layer_num = len(hidden_size_list)
        self.use_dropout = use_dropout
        self.weight_decay_lambda = weight_decay_lambda
        self.use_batchnorm = use_batchnorm
        self.params = {}

        # 重みの初期化
        self.__init_weight(weight_init_std)

        # レイヤの生成
        activation_layer = {'sigmoid': layers.Sigmoid, 'relu': layers.Relu}
        self.layers = OrderedDict()
        for idx in range(1, self.hidden_layer_num+1):
            self.layers['Affine' + str(idx)] = layers.Affine(self.params['W' + str(idx)], self.params['b' + str(idx)])
            if self.use_batchnorm:
                self.params['gamma' + str(idx)] = np.ones(hidden_size_list[idx-1])
                self.params['beta' + str(idx)] = np.zeros(hidden_size_list[idx-1])
                self.layers['BatchNorm' + str(idx)] = layers.BatchNormalization(self.params['gamma' + str(idx)], self.params['beta' + str(idx)])
                
            self.layers['Activation_function' + str(idx)] = activation_layer[activation]()
            
            if self.use_dropout:
                self.layers['Dropout' + str(idx)] = layers.Dropout(dropout_ratio)

        idx = self.hidden_layer_num + 1
        self.layers['Affine' + str(idx)] = layers.Affine(self.params['W' + str(idx)], self.params['b' + str(idx)])

        self.last_layer = layers.SoftmaxWithLoss()

    def __init_weight(self, weight_init_std):
        all_size_list = [self.input_size] + self.hidden_size_list + [self.output_size]
        for idx in range(1, len(all_size_list)):
            scale = weight_init_std
            if str(weight_init_std).lower() in ('relu', 'he'):
                scale = np.sqrt(2.0 / all_size_list[idx - 1])  # ReLUを使う場合に推奨される初期値
            elif str(weight_init_std).lower() in ('sigmoid', 'xavier'):
                scale = np.sqrt(1.0 / all_size_list[idx - 1])  # sigmoidを使う場合に推奨される初期値
            self.params['W' + str(idx)] = scale * np.random.randn(all_size_list[idx-1], all_size_list[idx])
            self.params['b' + str(idx)] = np.zeros(all_size_list[idx])

    def predict(self, x, train_flg=False):
        for key, layer in self.layers.items():
            if "Dropout" in key or "BatchNorm" in key:
                x = layer.forward(x, train_flg)
            else:
                x = layer.forward(x)

        return x

    def loss(self, x, d, train_flg=False):
        y = self.predict(x, train_flg)

        weight_decay = 0
        for idx in range(1, self.hidden_layer_num + 2):
            W = self.params['W' + str(idx)]
            weight_decay += 0.5 * self.weight_decay_lambda * np.sum(W**2)

        return self.last_layer.forward(y, d) + weight_decay

    def accuracy(self, X, D):
        Y = self.predict(X, train_flg=False)
        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

    def gradient(self, x, d):
        # forward
        self.loss(x, d, train_flg=True)

        # backward
        dout = 1
        dout = self.last_layer.backward(dout)

        layers = list(self.layers.values())
        layers.reverse()
        for layer in layers:
            dout = layer.backward(dout)

        # 設定
        grads = {}
        for idx in range(1, self.hidden_layer_num+2):
            grads['W' + str(idx)] = self.layers['Affine' + str(idx)].dW + self.weight_decay_lambda * self.params['W' + str(idx)]
            grads['b' + str(idx)] = self.layers['Affine' + str(idx)].db

            if self.use_batchnorm and idx != self.hidden_layer_num+1:
                grads['gamma' + str(idx)] = self.layers['BatchNorm' + str(idx)].dgamma
                grads['beta' + str(idx)] = self.layers['BatchNorm' + str(idx)].dbeta

        return grads

# バッチ正則化 layer
class BatchNormalization:
    '''
    gamma: スケール係数
    beta: オフセット
    momentum: 慣性
    running_mean: テスト時に使用する平均
    running_var: テスト時に使用する分散
    '''
    def __init__(self, gamma, beta, momentum=0.9, running_mean=None, running_var=None):
        self.gamma = gamma
        self.beta = beta
        self.momentum = momentum
        self.input_shape = None

        self.running_mean = running_mean
        self.running_var = running_var  
        
        # backward時に使用する中間データ
        self.batch_size = None
        self.xc = None
        self.std = None
        self.dgamma = None
        self.dbeta = None

    def forward(self, x, train_flg=True):
        if self.running_mean is None:
            N, D = x.shape
            self.running_mean = np.zeros(D)
            self.running_var = np.zeros(D)
                        
        if train_flg:
            mu = x.mean(axis=0) # 平均
            xc = x - mu # xをセンタリング
            var = np.mean(xc**2, axis=0) # 分散
            std = np.sqrt(var + 10e-7) # スケーリング
            xn = xc / std
            
            self.batch_size = x.shape[0]
            self.xc = xc
            self.xn = xn
            self.std = std
            self.running_mean = self.momentum * self.running_mean + (1-self.momentum) * mu # 平均値の加重平均
            self.running_var = self.momentum * self.running_var + (1-self.momentum) * var #分散値の加重平均
        else:
            xc = x - self.running_mean
            xn = xc / ((np.sqrt(self.running_var + 10e-7)))
            
        out = self.gamma * xn + self.beta 
        
        return out

    def backward(self, dout):
        dbeta = dout.sum(axis=0)
        dgamma = np.sum(self.xn * dout, axis=0)
        dxn = self.gamma * dout
        dxc = dxn / self.std
        dstd = -np.sum((dxn * self.xc) / (self.std * self.std), axis=0)
        dvar = 0.5 * dstd / self.std
        dxc += (2.0 / self.batch_size) * self.xc * dvar
        dmu = np.sum(dxc, axis=0)
        dx = dxc - dmu / self.batch_size
        
        self.dgamma = dgamma
        self.dbeta = dbeta

        return dx    
    
# データの読み込み
# (x_train, d_train), (x_test, d_test) = load_mnist(normalize=True, one_hot_label=True)
(x_train, d_train), (x_test, d_test) = load_mnist(normalize=True)

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

activations = ['sigmoid', 'relu']
weight_init_stds = [0.01, 'Xavier', 'He']
use_batchnorms = [False, True]

iters_num = 2000
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

plot_interval = 100
plot_idx = 0 

for k in range(len(activations)):
    for l in range(len(weight_init_stds)):
            for m in range(len(use_batchnorms)):
                network = MultiLayerNet(input_size=784, hidden_size_list=[40, 20], output_size=10, activation=activations[k], weight_init_std=weight_init_stds[l], use_batchnorm=use_batchnorms[m])

                train_loss_list = []
                accuracies_train = []
                accuracies_test = []
                lists = []
                plot_idx = plot_idx + 1

                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]

                    # 勾配
                    grad = network.gradient(x_batch, d_batch)

                    for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
                        network.params[key] -= learning_rate * grad[key]

                    loss = network.loss(x_batch, d_batch)
                    train_loss_list.append(loss)

                    if (i + 1) % plot_interval == 0:
                        accr_test = network.accuracy(x_test, d_test)
                        accuracies_test.append(accr_test)        
                        accr_train = network.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))


                lists = range(0, iters_num, plot_interval)

                plt.rcParams['figure.figsize'] = (12.0, 10.0)
                plt.subplot(4,3,plot_idx)
                plt.plot(lists, accuracies_train, label='training set')
                plt.plot(lists, accuracies_test,  label='test set')
                plt.legend(loc='lower right')
                plt.title(activations[k] + ', '  + str(weight_init_stds[l]) + ', バッチ正規化'  + str(use_batchnorms[m]) + ' ('  + str(np.round(accuracies_test[-1],2)) + ')')
                plt.xlabel('count')
                plt.ylabel('accuracy')
                plt.ylim(0, 1.0)
        
# グラフの表示
plt.tight_layout()
# plt.suptitle('活性化関数および重みの初期値を変更した場合の予測精度', fontsize = 16)
plt.show()

#Section2：学習率最適化手法
学習率の値が大きい場合、最適値にいつまでもたどり着かず発散してしまう。
学習率の値が小さい場合、発散することはないが、小さすぎると収束するまでに時間がかかったり、大域局所最適値に収束しづらくなってしまったりする。

• モメンタム
  $$ V_t = \mu V_{t-1}-\epsilon\nabla E $$
  $$ w^{(t+1)} = w^{(t)}+V_t $$
  誤差をパラメータで微分したものと学習率の積を減算した後、現在の重みに前回の重みを減算した値と慣性の積を加算する。

【モメンタムのメリット】
- 局所最適解にはならず、大域的最適解となる。
- 谷間に着いてから最も低い位置（最適値）に行くまでの時間が早い。

モメンタムの勾配

# 勾配
grad = network.gradient(x_batch, d_batch)
if i == 0:
    v = {}
for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
    if i == 0:
        v[key] = np.zeros_like(network.params[key])
    v[key] = momentum * v[key] - learning_rate * grad[key]
    network.params[key] += v[key]

    loss = network.loss(x_batch, d_batch)
    train_loss_list.append(loss)

• AdaGrad
  $$ h_0 = \theta $$
  $$ h_t = h_{t-1}+(\nabla E)^2 $$
  $$ w^{(t+1)} = w^{(t)}-\epsilon \frac{1}{\sqrt{h_t}+\theta}\nabla E $$
  誤差をパラメータで微分したものと再定義した学習率の積を減算する。

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

AdaGradの勾配

# 勾配
grad = network.gradient(x_batch, d_batch)
if i == 0:
    h = {}
for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
    if i == 0:
        h[key] = np.full_like(network.params[key], 1e-4)
    else:
        h[key] += np.square(grad[key])
    network.params[key] -= learning_rate * grad[key] / (np.sqrt(h[key]))

    loss = network.loss(x_batch, d_batch)
    train_loss_list.append(loss)    

• RMSProp
  $$ h_t = \alpha h_{t-1}+(1-\alpha)(\nabla E)^2 $$
  $$ w^{(t+1)} = w^{(t)}-\epsilon \frac{1}{\sqrt{h_t}+\theta}\nabla E $$
  誤差をパラメータで微分したものと再定義した学習率の積を減算する。

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

RMSPropの勾配

# 勾配
grad = network.gradient(x_batch, d_batch)
if i == 0:
    h = {}
for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
    if i == 0:
        h[key] = np.zeros_like(network.params[key])
    h[key] *= decay_rate
    h[key] += (1 - decay_rate) * np.square(grad[key])
    network.params[key] -= learning_rate * grad[key] / (np.sqrt(h[key]) + 1e-7)

    loss = network.loss(x_batch, d_batch)
    train_loss_list.appen

• Adam
  $$ V_t = \mu V_{t-1}-\epsilon\nabla E $$
  $$ h_t = \alpha h_{t-1}+(1-\alpha)(\nabla E)^2 $$
  $$ w^{(t+1)} = w^{(t)}-\epsilon \frac{V_t}{\sqrt{h_t}+\theta}\nabla E $$

【Adamのメリット】
- モメンタムの過去の勾配の指数関数的減衰平均、RMSPropの過去の勾配の２乗の指数関数的減衰平均というメリットを孕んだアルゴリズムである。

Adamの勾配

# 勾配
grad = network.gradient(x_batch, d_batch)
if i == 0:
    m = {}
    v = {}
learning_rate_t  = learning_rate * np.sqrt(1.0 - beta2 ** (i + 1)) / (1.0 - beta1 ** (i + 1))    
for key in ('W1', 'W2', 'W3', 'b1', 'b2', 'b3'):
    if i == 0:
        m[key] = np.zeros_like(network.params[key])
        v[key] = np.zeros_like(network.params[key])

    m[key] += (1 - beta1) * (grad[key] - m[key])
    v[key] += (1 - beta2) * (grad[key] ** 2 - v[key])            
    network.params[key] -= learning_rate_t * m[key] / (np.sqrt(v[key]) + 1e-7)                

    loss = network.loss(x_batch, d_batch)
    train_loss_list.append(loss)

#Section3：過学習
テスト誤差と訓練誤差とで学習曲線が乖離し、特定の訓練サンプルに対して特化して学習してしまうこと。
過学習を防ぐには、下記のような方法がある。

• L2ノルムを利用:Ridge推定量（縮小推定…パラメータを0に近づけるよう推定）
  $$ \sum_{i=1}^n(y_i-\beta_0-\sum_{j=1}^p\beta_jx_{ij})^2 + \lambda\sum_{j=1}^p\beta_j^2 $$
  $\lambda$：ハイパーパラメータ

  • $\lambda$がゼロであれば最小二乗法と同じ
  • $\lambda$を大きくすると$\beta_1,…,\beta_p$は0に近づく（$\beta_0$にはペナルティがつかないことに注意）
  • 交差検証などによって適切な値に決める

• L1ノルムを利用:Lasso（the Least absolute shrinkage and selection operator）推定量（スパース推定…いくつかのパラメータを正確に0に推定）
  $$ \sum_{i=1}^n(y_i-\beta_0-\sum_{j=1}^p\beta_jx_{ij})^2 + \lambda\sum_{j=1}^p|\beta_j| $$

  • パラメータのL1ノルムに比例するペナルティ
  • $\lambda$ を大きくすると多くのパラメータが0になる
    　$\rightarrow$変数選択によるスパースなモデルを生成

• ドロップアウト
  ランダムにノードを削除して学習させることにより、データ量を変化させずに異なるモデルを学習させていると解釈できる。

＜正則化なし（過学習の再現）＞
(optimizer.SGD(learning_rate=0.01))

＜L2正則化＞
(learning_rate=0.01)

＜L1正則化＞
(learning_rate=0.1)

＜ドロップアウト＞
(optimizer.SGD(learning_rate=0.01), weight_decay_lambda = 0.01)

(optimizer.Momentum(learning_rate=0.01, momentum=0.9), weight_decay_lambda = 0.01)

(optimizer.AdaGrad(learning_rate=0.01), weight_decay_lambda = 0.01)

(optimizer.Adam(learning_rate=0.01), weight_decay_lambda = 0.01)

＜ドロップアウト + L1正則化＞
(dropout_ratio = 0.1, weight_decay_lambda=0.005)

#Section4：畳み込みニューラルネットワークの概念

####畳み込みの演算概念

• バイアス
  

• （ゼロ）パディング
  

• ストライド
  

• チャンネル
  

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

全結合層のデメリット：画像の場合、縦、横、チャンネルの3次元データだが、1次元のデータとして処理される。すなわち、RGBの各チャンネル間の関連性が学習に反映されない。

#Section5：最新のCNN（2020年現在は最新とは呼べないが…）
####AlexNet

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

ChainerのAlexNetは、下記のようなコードになっている。

alex.py

import chainer
import chainer.functions as F
import chainer.links as L


class Alex(chainer.Chain):

    """Single-GPU AlexNet without partition toward the channel axis."""

    insize = 227

    def __init__(self):
        super(Alex, self).__init__()
        with self.init_scope():
            self.conv1 = L.Convolution2D(None, 96, 11, stride=4)
            self.conv2 = L.Convolution2D(None, 256, 5, pad=2)
            self.conv3 = L.Convolution2D(None, 384, 3, pad=1)
            self.conv4 = L.Convolution2D(None, 384, 3, pad=1)
            self.conv5 = L.Convolution2D(None, 256, 3, pad=1)
            self.fc6 = L.Linear(None, 4096)
            self.fc7 = L.Linear(None, 4096)
            self.fc8 = L.Linear(None, 1000)

    def __call__(self, x, t):
        h = F.max_pooling_2d(F.local_response_normalization(
            F.relu(self.conv1(x))), 3, stride=2)
        h = F.max_pooling_2d(F.local_response_normalization(
            F.relu(self.conv2(h))), 3, stride=2)
        h = F.relu(self.conv3(h))
        h = F.relu(self.conv4(h))
        h = F.max_pooling_2d(F.relu(self.conv5(h)), 3, stride=2)
        h = F.dropout(F.relu(self.fc6(h)))
        h = F.dropout(F.relu(self.fc7(h)))
        h = self.fc8(h)

        loss = F.softmax_cross_entropy(h, t)
        chainer.report({'loss': loss, 'accuracy': F.accuracy(h, t)}, self)
        return loss