🧠

【Python】ADALINE(フルバッチ勾配降下法)と学習の収束

2024/03/03に公開

scikit-learn

はじめに

この記事は株式会社インプレスの「Python機械学習プログラミング Pytorch&scilit-learn編」を読んで、私が学習したことをまとめています。リンク一覧はこちらから。
今回は２章３節のパーセプトロンの改良版アルゴリズムである「ADALINE(フルバッチ勾配降下法)と学習の収束」を読んで学んだことをまとめていきます。
なお、この記事の続きは【Python】ADALINE(確率的勾配降下法)と学習の収束にまとめていますので、よければご覧ください。

※本を読んでいてわかりづらいところが多々ありました(無定義な用語があったり、文字の意味が多義だったり、複数の名称が一つの定義を共有していたり)。数学畑の自分としてはそのまま読むと少し難しかったので、自分が読みやすくなるように定義を明文化することにしました。
　ただ、数学の分野と違って工学や情報学の分野は定義が明文化されているケースが少ないのか（私の探し方が不十分であることも十分考えられますが）、定義が書かれた参考文献を見つけることができませんでした。
　よってほぼゼロから自作した定義のため、やや正確性に欠ける拙い定義です。その都度修正をかけてより良いものにしたいと思いますが、まともな定義になるまで生温かい目で見てくださるとありがたいです。また、「この本を参考にするといいよ！」などアドバイスをいただけると、物凄く喜びます。

定義

データベースにおけるテーブルのことをデータセットと定義します。

あるデータセットにおいて、モデルのパラメータを適応的に調整する変数が並んだカラム(列)を特徴量(feature)といい、 $\boldsymbol{x}_{i}$ と書きます。ただし、 $i$ は特徴量の種類を表す自然数とします。また、カラムの変数がインデックス値などにより順序が与えられるとき、 $i$ 種類目の特徴量の中で $j$ 番目（ $j$ は自然数）の変数を $x_{i}^{(j)}$ と書きます。

学習したモデルにより予測される値を予測値といい、 $\hat{y}$ と書きます。特に特徴量の変数 $x_{i}^{(j)}$ に対応する予測値は $\hat{y}^{(j)}$ と書きます。予測ではなく $x_{i}^{(j)}$ に対応して観測によって与えられた値を正解値といい、 $y^{(j)}$ と書きます。正解値を含むカラムを目的変数(target)といいます。

特徴量 $x_1, x_2, \cdots, x_m$ が与えられたとき、次の式

\begin{equation*} w_1 x_1 + w_2 x_2 + \cdots + w_m x_m + b \ (w_1, w_2, \cdots, w_m, b \in \boldsymbol{R}) \end{equation*}

と書けるとき、この式を総入力(net input)といい、 $z$ で表します。また、 $w_1, w_2, \cdots, w_m$ を $x_1, x_2, \cdots, x_m$ それぞれに対応する重みといい、 $b$ をバイアスといいます。

つまり、行列表記なども用いると以下のように表すことができます。

\begin{alignat*}{2} z &= w_1 x_1 + w_2 x_2 + \cdots + w_m x_m + b \\ &= (w_1 \ w_2 \cdots \ w_m) \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_m \end{array} \right) + b \\ &= \boldsymbol{w}^{T} \boldsymbol{x} + b \end{alignat*}

ただし、

\begin{equation*} \boldsymbol{w} = \left( \begin{array}{c} w_1 \\ w_2 \\ \vdots \\ w_m \end{array} \right) , \boldsymbol{x} = \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_m \end{array} \right) \end{equation*}

とし、 $\boldsymbol{w}^{T}$ は $\boldsymbol{w}$ の転置行列を表します。

さて、ここまで定義してやっと本題のADALINEについて入っていこうと思います。
ADALINEは、学習過程で損失関数を最適化しますが、その最適化する損失関数 $L_1(\boldsymbol{w}, b)$ を以下のように定義します。

\begin{equation*} L_1(\boldsymbol{w}, b) = \cfrac{1}{n} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)}))^2 \end{equation*}

※参考書では $\cfrac{1}{2n}$ となっていますが、便宜上つけただけらしく、論理的には不都合なのでこの記事では分母の $2$ を外して議論します。
ただし $n$ は特徴量の変数の個数を表し、 $\sigma_{2}$ は線形活性化関数とします。

なおこの損失関数 $L_1(\boldsymbol{w}, b)$ は、平均二乗誤差(mean squared error:MSE)とも呼ばれています。

ADALINEのアルゴリズムでは、この損失関数が凸関数であることを利用して勾配降下法(gradient descent)という手法を用いて、重みとバイアスを最適化させます。

このとき、 $\boldsymbol{w}$ と $b$ を以下のように定義して更新していきます。

\begin{alignat*}{2} \boldsymbol{w} &:= \boldsymbol{w} + \Delta \boldsymbol{w} \\ b &:= b + \Delta b \end{alignat*}

ただし、 $\Delta \boldsymbol{w}, \Delta b$ をそれぞれ次のように定義します。

\begin{alignat*}{2} \Delta \boldsymbol{w} &= \eta \times \left( - \nabla_{w} L_1(\boldsymbol{w}, b) \right) \\ \Delta b &= \eta \times \left( - \nabla_{b} L_1(\boldsymbol{w}, b) \right) \end{alignat*}

ただし、 $\nabla_{w} L_1(\boldsymbol{w}, b)$ は $L_1(\boldsymbol{w}, b)$ の $w$ に関する偏導関数を表し、 $\nabla_{b} L_1(\boldsymbol{w}, b)$ は $L_1(\boldsymbol{w}, b)$ の $b$ に関する偏導関数を表します。

※参考書ではもう少しスマートな式の形で表されていますが、参考書で表記されている定義に合わせてここでは表現してあります。

$\nabla_{w} L_1(\boldsymbol{w}, b)$ については、重み $w_k \ (1 \le k \le m)$ ごとに偏微分します。計算は以下の通りです。
※参考書よりかなり細かく書いています。

\begin{alignat*}{2} \cfrac{\partial}{\partial w_k} L_1(\boldsymbol{w}, b) &= \cfrac{\partial}{\partial w_k} \cfrac{1}{n} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)}))^2 \\ &= \cfrac{1}{n} \cfrac{\partial}{\partial w_k} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)}))^2 \\ &= \cfrac{1}{n} \cfrac{\partial}{\partial w_k} \left( (y^{(1)} - \sigma_{2}(z^{(1)}))^2 + (y^{(2)} - \sigma_{2}(z^{(2)}))^2 + \cdots + (y^{(n)} - \sigma_{2}(z^{(n)}))^2 \right) \\ &= \cfrac{1}{n} \left( \cfrac{\partial}{\partial w_k}(y^{(1)} - \sigma_{2}(z^{(1)}))^2 + \cfrac{\partial}{\partial w_k}(y^{(2)} - \sigma_{2}(z^{(2)}))^2 + \cdots + \cfrac{\partial}{\partial w_k}(y^{(n)} - \sigma_{2}(z^{(n)}))^2 \right) \\ &= \cfrac{1}{n} \left( 2(y^{(1)} - \sigma_{2}(z^{(1)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(1)} - \sigma_{2}(z^{(1)}))} + 2(y^{(2)} - \sigma_{2}(z^{(2)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(2)} - \sigma_{2}(z^{(2)}))} + \cdots + 2(y^{(n)} - \sigma_{2}(z^{(n)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(n)} - \sigma_{2}(z^{(n)}))} \right) \\ \end{alignat*}

ここで下線部の式についてですが、普通はこのような形で計算は止めません。ただ、詳しく解説しようと思いあえて止めてあります。最初の下線部に注目してみると、次のように計算ができます。

\begin{alignat*}{2} \cfrac{\partial}{\partial w_k} (y^{(1)} - \sigma_{2}(z^{(1)})) &= \cfrac{\partial}{\partial w_k} y^{(1)} - \cfrac{\partial}{\partial w_k} \sigma_{2}(z^{(1)}) \\ &= 0 - \cfrac{\partial}{\partial w_k} (w_{1} x_{1}^{(1)} + w_{2} x_{2}^{(1)} + \cdots + \underline{w_{k} x_{k}^{(1)}} + \cdots + w_{n} x_{m}^{(1)} + b) \\ &= 0 - (0 + 0 + \cdots + \underline{x_{k}^{(1)}} + \cdots + 0 + 0) \\ &= - x_{k}^{(1)} \end{alignat*}

したがって、

\begin{alignat*}{2} \cfrac{\partial}{\partial w_k} L_1(\boldsymbol{w}, b) &= \cfrac{1}{n} \left( 2(y^{(1)} - \sigma_{2}(z^{(1)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(1)} - \sigma_{2}(z^{(1)}))} + 2(y^{(2)} - \sigma_{2}(z^{(2)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(2)} - \sigma_{2}(z^{(2)}))} + \cdots + 2(y^{(n)} - \sigma_{2}(z^{(n)})) \underline{\cfrac{\partial}{\partial w_k} (y^{(n)} - \sigma_{2}(z^{(n)}))} \right) \\ &= \cfrac{1}{n} \left( 2(y^{(1)} - \sigma_{2}(z^{(1)})) (- x_{k}^{(1)}) + 2(y^{(2)} - \sigma_{2}(z^{(2)})) (- x_{k}^{(2)}) + \cdots + 2(y^{(n)} - \sigma_{2}(z^{(n)})) (- x_{k}^{(n)}) \right) \\ &= \cfrac{2}{n} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)})) (- x_{k}^{(j)}) \end{alignat*}

となります。同様にして $\nabla_{b} L_1(\boldsymbol{w}, b)$ も $b$ で偏微分し、以下の結果を得ることができます。

\begin{alignat*}{2} \cfrac{\partial}{\partial b} L_1(\boldsymbol{w}, b) &= - \cfrac{2}{n} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)})) \end{alignat*}

以上より、重みの更新とバイアスの更新の式は以下のように表すことができます。
※のちの実装コードとできるだけ同じ形に表現しますので、是非比較してみてください。

\begin{alignat*}{2} \Delta \boldsymbol{w} &= \eta \times \left( - \nabla_{w} L_1(\boldsymbol{w}, b) \right) \\ &= \eta \times \left( \cfrac{2}{n} \sum_{j = 1}^{n} (\boldsymbol{x}^{(j)})^{T} (y^{(j)} - \sigma_{2}(z^{(j)})) \right) \\ &= \sum_{j = 1}^{n} \eta \times \cfrac{2 \times (\boldsymbol{x}^{(j)})^{T} (y^{(j)} - \sigma_{2}(z^{(j)}))}{n} \\ \Delta b &= \eta \times \left( - \nabla_{b} L_1(\boldsymbol{w}, b) \right) \\ &= \eta \times \left( \cfrac{2}{n} \sum_{j = 1}^{n} (y^{(j)} - \sigma_{2}(z^{(j)})) \right) \\ &= \sum_{j = 1}^{n} \eta \times \cfrac{2 \times (y^{(j)} - \sigma_{2}(z^{(j)}))}{n} \\ \end{alignat*}

※なおこのパラメータの更新は、訓練データごとにパラメータを小刻みに更新するのではなく、訓練データセットのすべての訓練データに基づいて計算されます。書籍ではこの勾配降下法を特にフルバッチ勾配降下法(full batch gradient descent)と呼んでいます。
以上を踏まえたうえで、ADALINEの実装から学習率の最適化までを見ていきたいと思います。

ライブラリのインポート

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

ADALINE(フルバッチ勾配降下法)の実装

ADAptive LInear NEuron分類器を実装していきいます。

class AdalineGD:
    """
    パラメータ
    ----------------------------------------------------
    eta          : float.学習率
    n_iter       : int.訓練データの訓練回数
    random_state : int.重みを初期化するための乱数シード
    ----------------------------------------------------

    データ属性
    ----------------------------------------------------
    w_      : 1次元配列.適合後の重み
    b_      : スカラー.適合後のバイアス
    losses_ : リスト.各エポックでの損失関数の値
    ----------------------------------------------------
    """

    def __init__(self, eta = 0.01, n_iter = 50, random_state = 1):
        self.eta          = eta             # 学習率の初期化
        self.n_iter       = n_iter          # 訓練回数の初期化
        self.random_state = random_state    # 乱数シードを設定


    def fit(self, X, y):
        """
        パラメータ
        -----------------------------------------------------------------------------------------------
        X : {配列のようなデータ構造}. shape = [n_examples(訓練データの個数), n_features(特徴量の個数)]
        y : {配列のようなデータ構造}. shape = [n_examples(訓練データの個数)]
        -----------------------------------------------------------------------------------------------
        """
        rgen         = np.random.RandomState(self.random_state)
        self.w_      = rgen.normal(loc   = 0.0,
                                   scale = 0.01,
                                   size  = X.shape[1]
                                  )
        self.b_      = np.float_(0.)
        self.losses_ = []

        # 訓練回数分まで訓練データを反復
        for i in range(self.n_iter):
            net_input = self.net_input(X)
            output    = self.activation(net_input)
            errors    = y - output                    # 誤差の計算

            # モデルのパラメータ(w_, b_)の更新を行う
            # 重み(w_)の更新（学習率 * 重みごとの損失関数の偏微分係数）
            # なお、パーセプトロンではこの更新段階で予測値を用いていた。
            self.w_ += self.eta * (2.0 * X.T.dot(errors) / X.shape[0])
            # バイアス(b_)の更新（学習率 * 損失関数のバイアスについての偏微分係数）
            self.b_ += self.eta * (2.0 * errors.mean())

            # 損失関数の計算（シグマがついていない形。シグマの意味はループによって形成される）
            loss = (errors ** 2).mean()

            # 損失値を格納
            self.losses_.append(loss)
        return self

    # 総入力を計算
    def net_input(self, X):
        return np.dot(X, self.w_) + self.b_

    # 線形活性化関数の出力を計算
    def activation(self, X):
        return X   # ただの恒等関数

    # 1ステップ後のクラスラベルを返す
    def predict(self, X):
        return np.where(self.activation(self.net_input(X)) >= 0.5, 1, 0)

IrisデータセットでADALINEモデルの訓練

次にデータセットの準備を行います。

s = 'https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data'
print('From URL:', s)

# IrisデータセットをDaraFrameオブジェクトに直接読み込む
df = pd.read_csv(s,
                 header   = None,
                 encoding = 'utf-8')
df.tail()

散布図を用いてデータを可視化する

がく片の長さと花びらの長さを抽出し、散布図を用いてデータを可視化していきます。

# 1-100行目の1, 3行目(がく片の長さ、花びらの長さ)の抽出　※今回扱う特徴量は2つ
X = df.iloc[0 : 100, [0, 2]].values

# 1-100行目の目的変数の抽出
y = df.iloc[0 : 100, 4].values
# Iris-setosaを0, Iris-versicolorを1に変換
y = np.where(y == 'Iris-setosa', 0, 1)


# 品種setosaのプロット(緑の●)
plt.scatter(x      = X[:50, 0],
            y      = X[:50, 1],
            color  = '#3F9877',
            marker = 'o',
            label  = 'Setosa')
# 品種Versicolorのプロット(青の■)
plt.scatter(x      = X[50 : 100, 0],
            y      = X[50 : 100, 1],
            color  = '#003F8E',
            marker = 's',
            label  = 'Versicolor')

# 軸のラベルの設定
plt.xlabel('Sepal length [cm]')
plt.ylabel('Petal length [cm]')
# 凡例の設定
plt.legend(loc = 'upper left')

# グラフを保存
plt.savefig('fig2-6.png')

plt.show()

ADALINEのアルゴリズムを用いて訓練する

2つの学習率を用いて、エポック数に対する損失関数をプロットします。

fig, ax = plt.subplots(nrows = 1,
                       ncols = 2,
                       figsize = (10, 4)
                       )

# 勾配降下法によるADALINEの学習(学習率0.1)
ada1 = AdalineGD(n_iter = 15, eta = 0.1).fit(X, y)

# エポック数と孫実関数の関係を表す折れ線グラフを描画する.なお、縦軸の損失関数は常用対数
ax[0].plot(range(1, len(ada1.losses_) + 1), np.log10(ada1.losses_), marker = 'o')
# 軸ラベルの設定
ax[0].set_xlabel('Epochs')
ax[0].set_ylabel('log(Mean squared error)')
# タイトルの設定
ax[0].set_title('Adaline - Learning rate 0.1')


# 勾配降下法によるADALINEの学習(0.0001)
ada2 = AdalineGD(n_iter = 15, eta = 0.0001).fit(X, y)

# エポック数と損失関数の関係を表す折れ線グラフのプロット
ax[1].plot(range(1, len(ada2.losses_) + 1), ada2.losses_, marker = 'o')
# 軸ラベルの設定
ax[1].set_xlabel('Epochs')
ax[1].set_ylabel('Mean squared error')
# タイトルの設定
ax[1].set_title('Adaline - Learning rate 0.0001')

# グラフを保存
plt.savefig('fig2-11.png')

plt.show()

不適切な学習率では、うまく損失関数を最適化（重みとバイアスの最適化）ができないことがわかります。
本来はこの学習率のようなハイパーパラメータの最適化を行うことが先決ですが、この後ある程度最適な値で学習の様子を描画していきます。その前にデータの標準化を行います。
標準化と勾配降下法の相性は良いことが知られており、学習率を探す際に効率よく探すことができるようになります（ただし今回は程よい学習率が判明している設定になっているので、ここでは直観的にその恩恵を感じることはできませんが）。

特徴量のスケーリング

# データのコピー
X_std = np.copy(X)

# 各列の標準化
X_std[:, 0] = (X[:, 0] - X[:, 0].mean()) / X[:, 0].std()
X_std[:, 1] = (X[:, 1] - X[:, 1].mean()) / X[:, 1].std()

決定境界を可視化するための関数の実装

2次元のデータセットの決定境界を可視化するための関数を実装します。

def plot_decision_regions(X, y, classifier, test_idx = None, resolution = 0.02):

    """マーカーとカラーマップの準備"""
    markers = ('o', 's', '^', 'v', '<')
    colors  = ('#3F9877',                  # ジェードグリーン
               '#003F8E',                  # インクブルー
               '#EA5506',                  # 赤橙
               'gray',
               'cyan'
              )
    cmap    = ListedColormap(colors[: len(np.unique(y))])

    """グリッドポイント（格子点）の生成"""    
    # x軸方向の最小値、最大値を定義
    x_min = X[:, 0].min() - 1
    x_max = X[:, 0].max() + 1
    # y軸方向の最小値、最大値を定義
    y_min = X[:, 1].min() - 1
    y_max = X[:, 1].max() + 1

    # 格子点の生成
    xx, yy = np.meshgrid(np.arange(x_min, x_max, resolution),
                         np.arange(y_min, y_max, resolution)
                        )
    # 確認用
    #print(xx1)
    #print(xx2)

    """各特徴量を1次元配列に変換(ravel())して予測を実行"""
    # つまり2つの特徴量から0と予測された格子点と1と予測された値が格子点ごとにlabに格納される
    lab = classifier.predict(np.array([xx.ravel(), yy.ravel()]).T)

    # 確認用
    #print(np.array([xx1.ravel(), xx2.ravel()]).T)
    #print(lab) # 格子点の並びで0,1が格納されているが、この時点では1次元配列なので変換が必要

    """予測結果の元のグリッドポイント(格子点)のデータサイズに変換"""
    lab = lab.reshape(xx.shape)

    """グリッドポイントの等高線のプロット"""
    plt.contourf(xx, yy, lab,
                 alpha = 0.3,                    # 透過度を指定
                 cmap  = cmap
                 )
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())

    """決定領域のプロット"""
    # クラスごとに訓練データをセット
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x         = X[y == cl, 0],
                    y         = X[y == cl, 1],
                    alpha     = 0.8,
                    c         = colors[idx],
                    marker    = markers[idx],
                    label     = f'Class {cl}',
                    edgecolor = 'black'
                    )
        
    """テストデータ点を目立たせる（点を〇で表示）"""
    if test_idx:    # ここはbool(test_idx)と同義。つまりTrueを返す
        # すべてのデータ点を描画
        X_test = X[test_idx, :]
        y_test = y[test_idx]

        plt.scatter(X_test[:, 0],
                    X_test[:, 1],
                    c         = 'none',
                    edgecolor = 'black',
                    alpha     = 1.0,
                    linewidth = 1,
                    marker    = 'o',
                    s         = 100,
                    label     = 'Test set'
                   )

最適な学習率を用いて再度学習し、決定境界を可視化

標準化の後、学習率η=0.5と、小さなエポック数にもとづいてADALINEを再び訓練し、収束することを確認します。

# 勾配降下法によるADALINEの学習(標準化後、学習率eta = 0.5)
ada_gd = AdalineGD(n_iter = 20, eta = 0.5)

# モデルの学習
ada_gd.fit(X_std, y)

# 決定領域のプロット
plot_decision_regions(X_std, y, classifier = ada_gd)

plt.title('Adaline - Gradient descent')
plt.xlabel('Sepal length [standardized]')
plt.ylabel('Petal length [standardized]')
plt.legend(loc = 'upper left')

plt.tight_layout()
plt.show()

# エポック数と損失値（MSE）の関係を表す折れ線グラフの描画
plt.plot(range(1, len(ada_gd.losses_) + 1), ada_gd.losses_, marker = 'o')

plt.xlabel('Epochs')
plt.ylabel('Mean squared error')
plt.tight_layout()

# グラフを保存
plt.savefig('fig2-14.png')

plt.show()