勾配降下法

勾配降下法(最急降下法)の具体的な計算についてまとまった記事が少なかったため、今回取り上げてみました。

勾配降下法とは

勾配降下法は、パラメータ(重みとバイアス)を少しずつ更新して勾配(傾き)が最小になるパラメータを探索する手法です。機械学習では、モデルが上手く学習しているかを計る手法の一つに、損失関数(予測した値と実際の値の差の関数)がありました。損失が少ないほど予測した値と実際の値の差は小さく、良いモデルだと言えます。ニューラルネットワーク等において、この損失関数に勾配降下法を適用することで、損失を最小化するようなパラメータを探し出すことができます。

例えば、以下の関数yを損失関数、xをパラメータとみなし、最初にランダムな点x = 3を選んだとしましょう。x = 3における勾配は4となります。このxをx = 2.5、x = 2、x = 1.5のように勾配が小さくなる方向へ進めると、やがてyが最小となるx = 1にたどり着きます。yが損失関数であることを思い出すと、損失が最小となるパラメータを探し当てたことになります。

具体例

早速ですが、勾配降下法の計算式は以下のようになります。「xから、勾配に学習率ηを掛けた分だけ差し引き、xを更新する」といった意味あいです。勾配がプラスの時はxがマイナス方向に更新され、勾配がマイナスの時はxがプラス方向に更新されます。学習率ηについては後で説明します。

\[
x := x – \eta\frac{df}{dx} \verb| ・・・①|
\]

以下の関数f(x)について、勾配降下法でf(x)が最小となるxを求めてみます。ちなみに、このグラフは変形して、\( f(x) = (x – 1) ^ 2 + 1 \)となるため、x = 1のときf(x)が最小値をとることがわかります。

\[
f(x) = x^2 – 2x + 2 \\
\]

xの初期値と学習率ηは、それぞれx0 = 3とη = 0.2の適当な値を選びました。本来左辺と右辺を同時更新していくので、添え字(x0、x1等の数字部分)はつけませんが、見やすさを優先してつけてみました。\( f(x)’ = 2x – 2 \)であることを考慮して、各値を①へ代入していくと、

\[
x_1 = x_0 – \eta\frac{df}{dx} = 3 – 0.2 \times (2 \times 3 – 2) = 2.2 \\
x_2 = x_1 – \eta\frac{df}{dx} = 2.2 – 0.2 \times (2 \times 2.2 – 2) = 1.72 \\
x_3 = x_2 – \eta\frac{df}{dx} = 1.72 – 0.2 \times (2 \times 1.72 – 2) = 1.43 \\
\]

このようにxを更新していくと、やがてx = 1へたどり着きます。手計算は大変なので、プログラムに頑張ってもらいましょう。後で使うライブラリもまとめてインポートしておきます。

# In
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Activation
# In
# 関数
def f(x):
    return (x - 1) ** 2 + 1

# 勾配
def grad(x):
    dx = 2 * x - 2 # fのxに関する微分
    return dx

# 勾配降下法
def gradient_descent(x, learning_rate=0.2, iteration=10): 
   
    # フレームを描画
    fig = plt.figure(figsize=(6, 6)) 
    ax = fig.add_subplot(111)
    plt.xticks(np.arange(-1, 5, step=1))
    plt.yticks(np.arange(0, 6, step=1))
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    
    # グラフを描画
    X = np.linspace(-1, 3)
    Y = f(X)
    ax.plot(X, Y, color='blue', linewidth=0.5)     
    
    # 初期値のx, yを描画
    ax.scatter(x, f(x), s=20, color='red', zorder=2) 
    
    # ループごとのx, yを描画
    for i in range(iteration):
        x =  x - (learning_rate * grad(x))
        ax.scatter(x, f(x), s=20, color='red', zorder=2)  
                
gradient_descent(3)

手計算とグラフのx座標が一致していることがわかります。10回程度探索を繰り返すとf(x)が最小となるx = 1へたどり着きました。上記では、1変数関数で探索しましたが、変数がいくつでも求めることができます。2変数関数だと以下のようになります。

# In
# 関数
def f(x, y):
    return x ** 2 + y ** 2

# 勾配と勾配の大きさ
def grad(x, y):
    dx = 2 * x # fのxに関する偏微分
    dy = 2 * y # fのyに関する偏微分
    norm = np.linalg.norm((dx, dy)) # 勾配の大きさ(ループ終了目的)
    return dx, dy, norm

# 勾配降下法
def gradient_descent(x, y, learning_rate=0.2, iteration=100):
    
    # フレームを描画
    fig = plt.figure(figsize=(10, 6)) 
    ax = fig.add_subplot(111, projection='3d') # 3DAxeSを追加
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    X, Y = np.mgrid[-10:10, -10:10] # 格子点の作成
    
    # グラフを描画
    Z = f(X, Y)
    ax.plot_wireframe(X, Y, Z, color='blue', zorder=1, linewidth=0.5) 
    ax.view_init(45, 45) # 視点の設定
    
    # 初期値のx, yを描画
    ax.scatter(x, y, f(x, y), s=20, color='red', zorder=2) 
    
    # ループごとのx, y, zを描画
    for i in range(1, iteration):
        x = x - learning_rate * grad(x, y)[0]
        y = y - learning_rate * grad(x, y)[1]
        ax.scatter(x, y, f(x, y), s=20, color='red', zorder=2)  

        # 勾配の大きさが微小になったらループ終了
        if grad(x, y)[2] < 1e-4: 
            break
            
gradient_descent(-10, -5)

学習率

学習率ηは、ランダムに選んだxから勾配(傾き)が最小になるxを探す際の、早さを表すハイパーパラメータ(定数)です。学習率ηが小さいと探索に時間がかかります。一方で、学習率ηが大きいと早く探索できますが、大きすぎると非効率が生じたり、逆に失敗したりします。

以下のグラフは、上記と同じ関数について学習率ηを変化させた時、xがどのようにふるまうかを表しています。

# In
# 関数
def f(x):
    return (x - 1) ** 2 + 1

# 勾配
def grad(x):
    dx = 2 * x - 2 # fのxに関する微分
    return dx

# 勾配降下法
def gradient_descent(x, learning_rates=[0.2, 0.5, 0.8, 1.05], iteration=10):
    
    # グラフフレームを描画
    fig = plt.figure(figsize=(10, 5))
    
    # 学習率をループ
    for learning_rate in learning_rates: 
        x_cord = [x] # xの初期値を格納(グラフ目的)
        iter = [i for i in range(iteration+1)] # イテレーション回数を格納(グラフ目的)
        
        # 学習率ごとに勾配降下法をループ
        for i in range(iteration): 
            x =  x - (learning_rate * grad(x))
            x_cord.append(x)
        
        # 学習率のループごとにグラフを描画
        plt.plot(iter, x_cord) 
        
        # 学習率のループごとにxを初期値(3)へ戻す
        x = x_cord[0] 
   
    plt.xticks(np.arange(0, iteration+1, step=1))
    plt.xlabel('Iteration')
    plt.ylabel('x')
    plt.legend(['η = 0.2', 'η = 0.5', 'η = 0.8', 'η = 1.05'])
    
gradient_descent(3)

損失関数への適用

「そんなまわりくどいことをしないで、損失関数を全てのパラメータに関して偏微分した勾配が0になるような点を見つければいいのでは?」というのはもっともな疑問です。1変数関数や2変数関数のように単純な関数の場合、簡単に損失関数の空間を計算(描画)できます。しかしながら、例えばニューラルネットワークの損失関数は何万ものパラメータを持っている(多次元空間でグラフがぐにゃぐにゃな)ケースがほとんどで、損失関数の空間を計算(描画)して、これが最小値を取るような(全てのパラメータに関する損失関数の偏微分が0になる)点を計算するには、かなりの計算リソースが必要です。勾配降下法を利用することで、こうした複雑な空間へ正面から立ち向かうのではなく、損失関数の値が減少するかどうかという1次元の線の問題へ転換することができるのです。

さて、損失関数の計算式は以下のようでした。iはi番目のデータを表しており、「1番目のデータからm番目のデータについて、予測した値と実際の値の差の平方を合計し、データ数で平均する」といった意味あいです。ここではxとyはパラメータでなく、それぞれ特徴量行列(定数)とターゲットベクトル(定数)であることを注意する必要があります。パラメータはwとbです。ちなみに、J(w, b)へ1/m(サンプル数の平均)ではなく1/2mをかけることで、合成関数の微分の過程で生じる2と相殺しています。J(w, b)の凹凸には影響しないので大きな影響はありませんが、勾配が緩やかになるため、学習率ηへこの影響を織り込む必要があります。

\[
J(w, b) = \sum^m_{i=1} \frac{1}{2m} (a(z^{(i)}) - y^{(i)}) ^2 \verb| ・・・②|\\
\]

a(z)は活性化関数で、特徴量行列とパラメータを掛け合わせた回帰モデルを想定しました。

\[
a(z) = w^Tx + b \verb| ・・・③|\\
\]

勾配降下法の計算式は以下のようでした。f(x)を損失関数J(w, b)に、xをパラメータwとbに置き換えています。ηはハイパーパラメータ(定数)ですから、パラメータwとbの偏微分さえわかれば、簡単に計算できます。

\[
w := w - \eta\frac{d}{dw}J(w, b) \verb| ・・・①(再掲)|\\
b := b - \eta\frac{d}{db}J(w, b) \verb| ・・・①(再掲)|
\]

パラメータwに関する偏微分

パラメータwに関して損失関数J(w, b)を偏微分します。⑤から⑥の式変形には、\( f(g(x)) = (a(z^{(i)}) - y^{(i)}) ^2 \)、\( g(x) = a(z^{(i)}) - y^{(i)} \)とし、合成関数の微分\( f(g(x))'=f'(g(x))g'(x) \)を利用しています。また、⑥から⑦の式変形には、後ろ半分のwに関する偏微分が\( (a(z^{(i)}) - y^{(i)})' = (w^Tx^{(i)} + b - y^{(i)})' = x^{(i)} \)となることを利用しています。

\[
\frac{d}{dw}J(w, b) = \frac{d}{dw} \sum^m_{i=1} \frac{1}{2m} (a(z^{(i)}) - y^{(i)}) ^2 \verb| ・・・④|\\
= \frac{1}{2m} \sum^m_{i=1} \frac{d}{dw} (a(z^{(i)}) - y^{(i)}) ^2 \verb| ・・・⑤|\\
= \frac{1}{2m} \sum^m_{i=1} 2 (a(z^{(i)}) - y^{(i)}) \frac{d}{dw} (a(z^{(i)}) - y^{(i)}) \verb| ・・・⑥|\\
= \frac{1}{m} \sum^m_{i=1} (a(z^{(i)}) - y^{(i)})x^{(i)} \verb| ・・・⑦|\\
\]

パラメータbに関する偏微分

パラメータbに関して損失関数J(w, b)を偏微分します。⑨から⑩の式変形には、上記と同様合成関数の微分を、⑩から⑪の式変形には、後ろ半分のbに関する偏微分が\( (a(z^{(i)}) - y^{(i)})' = (w^Tx^{(i)} + b - y^{(i)})' = 1 \)となることを利用しています。

\[
\frac{d}{db}J(w, b) = \frac{d}{db} \sum^m_{i=1} \frac{1}{2m} (a(z^{(i)}) - y^{(i)}) ^2 \verb| ・・・⑧|\\
= \frac{1}{2m} \sum^m_{i=1} \frac{d}{db} (a(z^{(i)}) - y^{(i)}) ^2 \verb| ・・・⑨|\\
= \frac{1}{2m} \sum^m_{i=1} 2 (a(z^{(i)}) - y^{(i)}) \frac{d}{db} (a(z^{(i)}) - y^{(i)}) \verb| ・・・⑩|\\
= \frac{1}{m} \sum^m_{i=1} (a(z^{(i)}) - y^{(i)}) \verb| ・・・⑪|\\
\]

これらの凄いところは、回帰モデルa(z)が複雑になっても機械的に損失関数J(w, b)の偏微分が求められることです。例えば、パラメータwの次数が増えても、⑦のxの次数を増やすだけで済みます。

最後に、⑦と⑪を①へ代入し、勾配降下法により損失関数J(w, b)が最小となるようなパラメータwとbを求めることができます。

実装

ボストン市の住宅価格のデータセットを利用し、損失及びパラメータwとbを求めたいと思います。まずはデータセットをインポートし、トレーニングデータとテストデータに分離後、スケーリングします。

# In
boston = datasets.load_boston()
X = boston.data
y = boston.target

# スケーリング前にトレーニングデータとテストデータを分離
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=123)

# スケーリング
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

損失の推移を観察するため、ニューラルネットワークモデルを用意しました。3層からなる全結合の単純な構成で、損失関数の評価には②で紹介した平均二乗誤差(MSE)を利用します。

# In
# 単純なニューラルネットワークモデルを作成
model = Sequential()
model.add(Dense(128, activation='relu'))
model.add(Dense(64, activation='relu'))
model.add(Dense(1))
model.compile(optimizer='adam',loss='mse')

トレーニングデータでモデルを学習させます。トレーニングデータを学習させる回数だけ、損失が計算されます。

# In
#トレーニングデータでモデルを学習
model.fit(x=X_train, 
          y=y_train, 
          epochs=100, # 学習回数
          validation_data=(X_test, y_test), # 今回はテストデータを検証データとする
          verbose=1)
# In
model.summary()

イテレーション(エポック)ごとの損失は以下のようになりました。トレーニングデータも、検証データも損失が減り続け、やがて収束していることから、過学習は起きていません。※今回はテストデータを検証データとしたことを注意下さい。

# In
# イテレーション(エポック)ごとの損失
losses = pd.DataFrame(model.history.history)
losses.plot(figsize=(6, 4))
plt.xlabel('Iteration')
plt.ylabel('Loss')

更新されたパラメータwとbを見ることができます。3層分のパラメータ約1万個が格納されているため、一部のみ抜粋します。

# In
for i in range(len(model.layers)):
    w = model.layers[i].get_weights()[0].T # 転置により活性化関数の形と合わせた(shape = (レイヤLのニューロン数, レイヤL-1のニューロン数))
    b = model.layers[i].get_weights()[1]
    print('=========================================================================')
    print('Layer {}: w (shape={})\n{}'.format(i+1, w.shape, w))
    print('-------------------------------------------------------------------------')
    print('Layer {}: b (shape={})\n{}'.format(i+1, b.shape, b))

このように1イテレーション(エポック)ごとに、現在のパラメータwとbで全トレーニングデータあるいは検証データの損失を計算・合計し、勾配降下法によりパラメータwとbを更新します。損失が減らなくなるまで繰り返すことで、損失を最小化するようなパラメータwとbを探し出すことができるのです。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です