勾配降下法(最急降下法)の具体的な計算についてまとまった記事が少なかったため、今回取り上げてみました。
勾配降下法とは
勾配降下法は、パラメータ(重みとバイアス)を少しずつ更新して勾配(傾き)が最小になるパラメータを探索する手法です。機械学習では、モデルが上手く学習しているかを計る手法の一つに、損失関数(予測した値と実際の値の差の関数)がありました。損失が少ないほど予測した値と実際の値の差は小さく、良いモデルだと言えます。ニューラルネットワーク等において、この損失関数に勾配降下法を適用することで、損失を最小化するようなパラメータを探し出すことができます。
例えば、以下の関数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を探し出すことができるのです。
勾配降下法は、パラメータ(重みとバイアス)を少しずつ更新して勾配(傾き)が最小になるパラメータを探索する手法です。機械学習では、モデルが上手く学習しているかを計る手法の一つに、損失関数(予測した値と実際の値の差の関数)がありました。損失が少ないほど予測した値と実際の値の差は小さく、良いモデルだと言えます。ニューラルネットワーク等において、この損失関数に勾配降下法を適用することで、損失を最小化するようなパラメータを探し出すことができます。
例えば、以下の関数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を探し出すことができるのです。
早速ですが、勾配降下法の計算式は以下のようになります。「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を探し出すことができるのです。
学習率ηは、ランダムに選んだ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を探し出すことができるのです。
「そんなまわりくどいことをしないで、損失関数を全てのパラメータに関して偏微分した勾配が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を探し出すことができるのです。
ボストン市の住宅価格のデータセットを利用し、損失及びパラメータ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を探し出すことができるのです。