Pythonで実装するニューラルネットワークを完成させよう：ニューラルネットワーク入門

ニューラルネットワークをNumPy（線形代数）なしのPythonでフルスクラッチ実装する連載（基礎編）の前／中／後編の後編。いよいよ完成。最適化の処理をPythonコードから理解しよう。

[一色政彦，デジタルアドバンテージ]

連載目次

　本稿は、ニューラルネットワーク（以下、ニューラルネット）の仕組みや挙動を、数学理論からではなくPythonコードから学ぶことを目標とした連載（基礎編）の最後となる第3回です。

　前々回の第1回では、「ニューラルネットの訓練（学習）処理を実現するために必要なこと」として、

• ステップ（1）順伝播： forward_prop()関数として実装（前々回）
• ステップ（2）逆伝播： back_prop()関数として実装（前回）
• ステップ（3）パラメーター（重みとバイアス）の更新： update_params()関数として実装（今回）。これによりモデルが最適化される

という3大ステップを示しました。前回の第2回で、このうちの「ステップ（2）逆伝播」までの実装が完了しています。

　今回はその続きとして、「ステップ（3）パラメーターの更新と、モデルの最適化」までを実装して、ニューラルネットの実装を完了させます。最後に、それを使って簡単な回帰問題を解いてみます。

　ここからの内容は簡単です。その分、作業的なコーディング部分が多くなってしまいますが、最後の完成に至るまでを楽しんでコーディングしていきましょう。

　※なお本稿は、第1回〜第2回とセットの内容なので、図や掲載コード（「リスト＜数字＞」と表記）などの番号は前回からの継続となっています。

ステップ3. パラメーター（重みとバイアス）更新の実装

　それでは、今回も何も見ずにゼロからスクラッチでコードを書くという想定で進めていきます。

　実装を始める前に、まずはもう一度、訓練（学習、最適化）処理全体の実装から振り返っておきましょう（前々回のリスト1）。

# 訓練処理
y_pred, cached_outs, cached_sums = forward_prop(cache_mode=True)  # （1）
grads_w, grads_b = back_prop(y_true, cached_outs, cached_sums)  # （2）
weights, biases = update_params(grads_w, grads_b, LEARNING_RATE)  # （3）

リスト1（抜粋して再掲）　訓練（学習）処理全体の実装

　前回は、back_prop()関数から戻り値としてgrads_w, grads_bという2つの勾配（gradients）情報を取得できるようにしました。今回はその勾配情報を使って、パラメーター（重みとバイアス）を更新していきます。

　パラメーターを更新する目的は、もちろん分かっていると思いますが、ニューラルネットのモデルを最適化することです。以下に、最適化の考え方を簡単にまとめておきます。

　※なお、本稿で説明するのは最も基礎的な勾配降下法（Gradient Descent）です。後述するSGD（確率的勾配降下法）もその一種で、他にはRMSpropやAdamなどより応用的な手法があります。SGD以外の場合は、重みパラメーターの更新方法も少し変わってきます。

　図14は、教材などでよく見る「最適化の参考イメージ」で、1つの重みパラメーターしかない場合の損失関数のグラフです。このようなイメージで、重みパラメーターの数値を調整（＝更新）しながら、坂の一番下まで進めていきます（※ちなみに、この図は2次元のグラフですが、2つの重みパラメーターがある場合の、3次元のグラフも教材などでよく見ます。3次元の場合は、谷の底に着くまで進めるイメージになります）。

図14　最適化の参考イメージ

　更新時に、どちらの方向に、どれくらい進むかを表す数値が、前回で計算済みの「勾配」となります。坂の下に進むには、現在の「各重み」から「その勾配」を引き算した新しい「各重み」に更新します（※例えば、勾配が−の数値の場合は、坂道は右下りになります。坂の下に進むには今の数値に＋しなければなりません。今の数値から−の勾配値を引く（−する）と＋になるので、坂道を右下に進めます。勾配が＋の数値の場合は、坂道を左下に進めます）。

　ただし1回の更新（イテレーション）で一気に進むと、曲線の中を左右に行ったり来たりしてなかなか収束しない可能性があります。そこで1回（＝オンライン学習なら1件のデータ、ミニバッチ学習ならバッチサイズ分のデータ）で進む距離が適度になるように、「学習率」（learning rate、η：イータ）を「勾配」に掛け算することで、進む大きさをスケーリングして調整します。

1つのパラメーターの更新

　1つの重みパラメーターの更新をPythonコードで書くと、リスト23のようになります。簡単ですね。バイアスの場合も、wをbに変えるだけの同じ式です。

# 取りあえず仮で、変数を定義して、コードが実行できるようにしておく
w_ij = 0.0  # 各重み
b_j = 0.0  # バイアス
grad_w_ij = 0.2  # 各重みの勾配
grad_b_j = 0.2  # バイアスの勾配
LEARNING_RATE = 0.1  # 学習率（lr）
lr = LEARNING_RATE
# ---ここまでは仮の実装。ここからが必要な実装---

w_ij = w_ij - lr * grad_w_ij  # 重みパラメーターの更新

b_j = b_j - lr * grad_b_j  # バイアスパラメーターの更新

リスト23　重みパラメーターの更新（SGDの場合）

　このコードは図15の計算式を表現したものです。

図15　重みパラメーター更新の計算式

　以上が分かれば、パラメーター（重みとバイアス）を更新するPython関数は実装できます。

パラメーター更新の処理全体の実装

　またまた順伝播の実装と同じ説明内容になりますが、ニューラルネットは、層があり、その中に複数のノードが存在するという構造ですので、

• 各層を1つずつ処理するforループと
  • 層の中のノードを1つずつ処理するforループの2段階構造が必要で
    • その中に「1つのパラメーターの更新」を記述

すればよいわけです。

　この考えに沿って、パラメーター更新の処理全体を行うupdate_params()関数を実装してみたのがリスト24です。2つのforループと、既に説明済みの「パラメーター更新」の実装（特に太字で示した4行）に注目してください。それら以外のコードは、更新した新しい「重み」と「バイアス」を多次元リストにまとめるためのこまごました処理なので、読み飛ばしても構いません。

def update_params(layers, weights, biases, grads_w, grads_b, lr=0.1):
    """
    パラメーター（重みとバイアス）を更新する関数。
    - 引数：
    (layers, weights, biases)： モデルを指定する。
    grads_w： 重みの勾配。
    grads_b： バイアスの勾配。
    lr： 学習率（learning rate）。最適化を進める量を調整する。
    - 戻り値：
    新しい重みとバイアスを返す。
    """

    # ネットワーク全体で勾配を保持するためのリスト
    new_weights = [] # 重み
    new_biases = [] # バイアス

    SKIP_INPUT_LAYER = 1
    for layer_i, layer in enumerate(layers):  # 各層を処理
        if layer_i == 0:
            continue  # 入力層はスキップ

        # 層ごとで勾配を保持するためのリスト
        layer_w = []
        layer_b = []

        for node_i in range(layer):  # 層の中の各ノードを処理
            b = biases[layer_i - SKIP_INPUT_LAYER][node_i]
            grad_b = grads_b[layer_i - SKIP_INPUT_LAYER][node_i]
            b = b - lr * grad_b  # バイアスパラメーターの更新
            layer_b.append(b)

            node_weights = weights[layer_i - SKIP_INPUT_LAYER][node_i]
            node_w = []
            for each_w_i, w in enumerate(node_weights):
                grad_w = grads_w[layer_i - SKIP_INPUT_LAYER][node_i][each_w_i]
                w = w - lr * grad_w  # 重みパラメーターの更新
                node_w.append(w)
            layer_w.append(node_w)

        new_weights.append(layer_w)
        new_biases.append(layer_b)
    
    return (new_weights, new_biases)

リスト24　パラメーター更新の処理全体の実装（SGDの場合）

　注意点は特にありません。コード中のコメントを参考にしてください。

　以上でupdate_params()関数が完成したので、試しに実行してみましょう。

パラメーター更新の実行例

　リスト25のようなコードを書けば、順伝播から逆伝播、パラメーター更新までを続けて実行できます。

layers = [2, 2, 2]
weights = [
    [[0.15, 0.2], [0.25, 0.3]],
    [[0.4, 0.45], [0.5,0.55]]
]
biases = [[0.35, 0.35], [0.6, 0.6]]
model = (layers, weights, biases)

# 元の重み
print(f'old-weights={weights}')
print(f'old-biases={biases}' )
# old-weights=[[[0.15, 0.2], [0.25, 0.3]], [[0.4, 0.45], [0.5, 0.55]]]
# old-biases=[[0.35, 0.35], [0.6, 0.6]]

# （1）順伝播の実行例
x = [0.05, 0.1]
y_pred, cached_outs, cached_sums = forward_prop(*model, x, cache_mode=True)

# （2）逆伝播の実行例
y_true = [0.01, 0.99]
grads_w, grads_b = back_prop(*model, y_true, cached_outs, cached_sums)
print(f'grads_w={grads_w}')
print(f'grads_b={grads_b}')
# grads_w=[[[0.006706025259285303, 0.013412050518570607], [0.007487461943833829, 0.014974923887667657]], [[0.6501681244277691, 0.6541291517796395], [0.13937181955411934, 0.1402209162240302]]]
# grads_b=[[0.13412050518570606, 0.14974923887667657], [1.09590596705977, 0.23492140409646534]]

# （3）パラメーター更新の実行例
LEARNING_RATE = 0.1 # 学習率（lr）
weights, biases = update_params(*model, grads_w, grads_b, lr=LEARNING_RATE)

# 更新後の新しい重み
print(f'new-weights={weights}')
print(f'new-biases={biases}')
# new-weights=[[[0.14932939747407145, 0.19865879494814295], [0.2492512538056166, 0.2985025076112332]], [[0.3349831875572231, 0.3845870848220361], [0.48606281804458806, 0.5359779083775971]]]
# new-biases=[[0.3365879494814294, 0.33502507611233234], [0.490409403294023, 0.5765078595903534]]

# モデルの最適化
model = (layers, weights, biases)

リスト25　パラメーター更新の実行例

　以上で、「ステップ（1）順伝播」「ステップ（2）逆伝播」「ステップ（3）パラメーター（重みとバイアス）の更新」を担う3つの関数の実装が完了しました。あとは、これら3つの関数を呼び出す最適化処理を実装して完成です。

3つのステップを呼び出す最適化処理の実装

最適化処理：学習方法と勾配降下法

　最適化処理を行う代表的な学習方法には幾つかの種類があります（参考：「ディープラーニング最速入門」）。実装を始める前に、代表的な学習方法を簡単にまとめておきます。

• オンライン学習（Online training）： データ1件ずつ訓練していくこと
• ミニバッチ学習（Mini-batch training）： 小さなまとまりのデータごとに訓練していくこと
• バッチ学習（Batch training）： データ全件で訓練していくこと

　このうち、オンライン学習やミニバッチ学習では、データをランダムにシャッフルすること（＝統計学の確率論におけるランダムサンプリング、無作為抽出をすること）で、（標本／サンプルである）訓練ごとのデータの分布が、（母集団である）データ全体の縮図になるようにします（※なお、バッチ学習のデータをシャッフルしても、データ全件を使うので無意味です。訓練時ではなく評価時は、再現性を保つためにもシャッフルしないのが一般的です）。このため、オンライン学習やミニバッチ学習の勾配降下法は、「確率的」勾配降下法（SGD）と呼ばれます。

　学習方法ごとに、勾配降下法をまとめると以下のようになります。

• SGD（Stochastic Gradient Descent）： オンライン学習
• ミニバッチSGD（Mini-batch SGD）： ミニバッチ学習。単にミニバッチ勾配降下法（Mini-batch Gradient Descent）とも呼ぶ
• 最急降下法（Steepest Descent）： バッチ学習。バッチ勾配降下法（Batch Gradient Descent）とも呼ぶ

　最適化の実装は、バッチサイズでデータを処理しない分、オンライン学習のSGDが一番簡単で、説明する上でも分かりやすい内容になると思います。しかし今回は、コード内容も簡単なので少し難易度を上げて、あえて全ての学習方法に対応できる実装コードにしてみます。

最適化の処理全体の実装

　訓練処理では、エポック（＝全データ分で1回の訓練）があり、その中にイテレーション（＝バッチサイズごとでのパラメーターの更新）が存在するという構造ですので、

• エポックを1回ずつ処理するforループと
  • その中にデータを1件ずつ処理するforループの2段階構造を用意し
    • その中に「ステップ（1）順伝播」「ステップ（2）逆伝播」と
    • イテレーションごとに「ステップ（3）パラメーターの更新」を記述

するようにします（※あくまで筆者による実装方針の例です）。

　この考えに沿って訓練（最適化）処理を実装しますが、階層が深くなる上にコードの行数が少し長いので、説明の都合上、上の箇条書きの前半2行をリスト26（train()親関数）、後半2行をリスト27（optimize()子関数）、という親子関係の2つの関数に分けて記述します。※1つの関数として実装した方がシンプルになって見通しもよくなるので、本来であればそうした方がよいと思います。

　まずリスト26の訓練処理では、2つのforループと、「訓練データのインデックスをランダムにシャッフル」している部分、リスト27で実装するoptimize()関数を呼び出して戻り値として返された損失値（loss）を蓄積（accumulate）している部分（acm_loss変数）（特に太字で示した7行）に注目してください。※シャッフルはデータ内容自体ではなく、データのインデックスだけをシャッフルしています。

import random

# 取りあえず仮で、空の関数を定義して、コードが実行できるようにしておく
def optimize(model, x, y, data_i, last_i, batch_i, batch_size, acm_g, lr=0.1):
    " モデルを最適化する関数（子関数）。"
    loss = 0.1
    return model, loss, batch_i, acm_g

# ---ここまでは仮の実装。ここからが必要な実装---

def train(model, x, y, batch_size=32, epochs=10, lr=0.1, verbose=10):
    """
    モデルの訓練を行う関数（親関数）。
    - 引数：
    model： モデルをタプル「(layers, weights, biases)」で指定する。
    x： 訓練データ（各データが行、各特徴量が列の、2次元リスト値）。
    y： 訓練ラベル（各データが行、各正解値が列の、2次元リスト値）。
    batch_size： バッチサイズ。何件のデータをまとめて処理するか。
    epochs： エポック数。全データ分で何回、訓練するか。
    lr： 学習率（learning rate）。最適化を進める量を調整する。
    verbose： 訓練状況を何エポックおきに出力するか。
    - 戻り値：
    損失値の履歴を返す。これを使って損失値の推移グラフが描ける。
    """
    loss_history = []  # 損失値の履歴

    data_size = len(y)  # 訓練データ数
    data_indexes = range(data_size)  # 訓練データのインデックス

    # 各エポックを処理
    for epoch_i in range(1, epochs + 1):  # 経過表示用に1スタート

        acm_loss = 0  # 損失値を蓄積（accumulate）していく

        # 訓練データのインデックスをシャッフル（ランダムサンプリング）
        random_indexes = random.sample(data_indexes, data_size)
        last_i = random_indexes[-1]  # 最後の訓練データのインデックス

        # 親関数で管理すべき変数
        acm_g = (None, None)  # 重み／バイアスの勾配を蓄積していくため
        batch_i = 0  # バッチ番号をインクリメントしていくため

        # 訓練データを1件1件処理していく
        for data_i in random_indexes:

            # 親子に分割したうちの子関数を呼び出す
            model, loss, batch_i, acm_g = optimize(
                model, x, y, data_i, last_i, batch_i, batch_size, acm_g, lr)

            acm_loss += loss  # 損失値を蓄積

        # エポックごとに損失値を計算。今回の実装では「平均」する
        layers = model[0]  # レイヤー構造
        out_count = layers[-1]  # 出力層のノード数
        # 「訓練データ数（イテレーション数×バッチサイズ）×出力ノード数」で平均
        epoch_loss = acm_loss / (data_size * out_count)

        # 訓練状況を出力
        if verbose != 0 and \
            (epoch_i % verbose == 0 or epoch_i == 1 or epoch_i == EPOCHS):
            print(f'[Epoch {epoch_i}/{EPOCHS}] train_loss: {epoch_loss}')

        loss_history.append(epoch_loss)  # 損失値の履歴として保存

    return model, loss_history


# サンプル実行用の仮のモデルとデータ
layers = [2, 2, 2]
weights = [
    [[0.15, 0.2], [0.25, 0.3]],
    [[0.4, 0.45], [0.5,0.55]]
]
biases = [[0.35, 0.35], [0.6, 0.6]]
model = (layers, weights, biases)
x = [[0.05, 0.1]]
y = [[0.01, 0.99]]

# モデルを訓練する
BATCH_SIZE = 2  # バッチサイズ
EPOCHS = 1  # エポック数
LEARNING_RATE = 0.02 # 学習率（lr）
model, loss_history = train(model, x, y, BATCH_SIZE, EPOCHS, LEARNING_RATE)
# 出力例：
# [Epoch 1/1] train_loss: 0.05

リスト26　訓練（最適化）処理全体を担う関数の実装

　コード中にもコメントを入れていますが、気を付けてほしいポイントを以下でも触れておきます。

　リスト26では、各データごとに蓄積した損失値（＝各エポック内の損失値の合計）をデータ数（data_size、出力ノード数のout_countは後述）で割って平均しています（厳密には、イテレーション、つまりバッチサイズごとの平均損失値を計算した後で、エポックごとに平均損失値を計算する方がよいです）。

　前回の実装では、損失関数に二乗和誤差（SSE：Sum of Squared Error）を採用しました。よって、合計の1/2（1/2 × Σ）という損失値を計算すべきです（1/2の部分は、損失関数やその偏導関数の実装の中で計算済みです）。しかし、リスト26でデータ数（上記の通り、厳密にはバッチサイズ）で平均している段階で、実質的には「1/2した平均二乗誤差（MSE：Mean Squared Error）」を使っていることに相当します（平均の1/2で1/2n × Σになっているので、完全なMSEでもないです）。

　このように損失関数を和（合計）から平均に変更したので、各勾配もバッチサイズごとに蓄積（後掲のリスト27のacm_gw変数やacm_gb変数）した後で、それをバッチサイズ（現在のバッチ数を意味するbatch_i）で割って平均する必要があります。

　なお、この「平均」の部分は、単に合計（Σ）や、合計の1/2（1/2 × Σ）、平均（1/n × Σ）、平均の1/2（1/2n × Σ）、その他の正規化など、のいずれのスケール調整（スケーリング）を行っても、最終的には学習率によってスケール調整されることになるので、結果は基本的に変わりません。しかし合計だけの場合、バッチサイズを変更するたびに学習率でスケールを調整しなければならなくなります。平均ならスケールは変わらないので調整しなくてよいです。そのため本稿では「平均」することにしました。

　実際にどうするかは、基本的に実装の目的や損失関数の計算式に依存します。分類問題では、交差エントロピー誤差などの損失関数に従い、基本的に平均します。また回帰問題では、平均二乗誤差などの損失関数に従い、基本的に平均します。