はじめに

以前、「機械学習の分類」で取り上げたアルゴリズムについて、その理論とpythonでの実装、scikit-learnを使った分析についてステップバイステップで学習していく。個人の学習用として書いてるので間違いなんかは大目に見て欲しいと思います。

今回は、損失関数のパラメータ最適化の手法である勾配降下法について。数学的に解が見つかる場合はよいが、世の中そういう場合の方が少ない。最適値を探索する手法としていろんな勾配方があるが、いくつかをコードを交えてまとめてみたい。例によって参考にしたサイトは以下です。ありがとうございます。

勾配法について

Wikipediaによると、勾配法とは

勾配法（こうばいほう、Gradient method）は、最適化問題において、関数の勾配に関する情報を解の探索に用いるアルゴリズムの総称。

と書いてある。単回帰では、$$ y = Ax+B $$という関数で近似させるために、損失関数$E$を残差の二乗和 $$ E = \sum_{i=0}^{N}(y_i - (Ax_i -B))^2 $$として、$E$を最小化するような$A$と$B$を求めた。この例では数学的に$A$と$B$を求めることができるのだが、損失関数がもっと複雑な場合であっても、最小となる（であろう）点を探すアプローチが勾配降下法である。

お題

例によってscikit-learnのdiabetes(糖尿病データ)を使います。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import datasets

diabetes = datasets.load_diabetes()

df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)

x = df['bmi'].values
y = diabetes.target
plt.scatter(x, y)

数学的に解いたときの傾きと切片はそれぞれ

傾きA:  949.43526038395
切片B:  152.1334841628967

でした。

最急勾配降下法

最急降下法(Steepest Gradient Descent)は、損失関数上のある点での勾配ベクトルと逆方向に少しずつ進んでいくというアイデアです。傾きを$\nabla E$とし、$$ E_{t+1}=E_{t} - \eta\nabla E_{t} $$を計算していきます。なお、$\eta$は学習率と呼ばれ、どの程度ゆっくり進んでいくかというものです。これが大きいと解が収束せずに発散してしまいます。

損失関数$$ E = \sum_{i=0}^{N}(y_i - (Ax_i -B))^2 $$について、$A$と$B$を求めるためには、$E$を$A$と$B$それぞれについて偏微分します。($\sum$の添え字は省略)

\frac{\partial{E}}{\partial{A}}=\sum(-2x_iy_i+2Ax_i^2+2Bx_i)
 \\
=-2\sum(Ax_i+B-y_i)x_i \\
\frac{\partial{E}}{\partial{B}}=\sum(-2x_iy_i+2Ax_i+2B)
 \\
=-2\sum(Ax_i+B-y_i) \\

となる。初期値を適当に決め、上の式を使って$A$と$B$を勾配方向に変化させていき、ほどよい段階で計算を打ち切ればよいことになる。

初期値の決め方、学習率の決め方、計算を打ち切る条件についてはさらなる考察が必要だが概略はこういうことだ。

pythonでの実装

最急降下法を実現するSteepestGradientDescentクラスを作成した。メソッドはなんとなくscikit-learnにあわせてみた。学習率は適当で、誤差が十分小さくなったらではなく回数で打ち切りとした。

class SteepestGradientDescent:
  def __init__(self, eta=0.001, n_iter=2000):
    self.eta = eta
    self.n_iter = n_iter
    self.grad = np.zeros((2,))
    self.loss = np.array([])

  def fit(self, X, Y, w0):
    self.w = w0
    for _ in range(self.n_iter):
      self.loss = np.append(self.loss, np.sum((Y-self.w[0]-self.w[1]*X)**2))
      self.grad[0] = np.sum(-2*(Y-self.w[0]-self.w[1]*X))
      self.grad[1] = np.sum(-2*X*(Y-self.w[0]-self.w[1]*X))

      self.w -= self.eta * self.grad

  def predict(self, x):
    return (self.w[0] + self.w[1]*x)

  @property
  def coef_(self):
    return self.w[1]

  @property
  def intercept_(self):
    return self.w[0]

  @property
  def loss_(self):
    return self.loss

fitメソッドの中で全データに対して勾配を算出し、係数を更新している。損失関数の値も同時に保存して、学習がどういう風に進んでいるかわかるようにした。

このクラスを使って糖尿病データに近似させ、損失の変化をグラフにプロットしてみる。

w0 = np.array([1.,1.])

model = SteepestGradientDescent()
model.fit(x, y, w0)

print("A: ", model.coef_)
print("B: ", model.intercept_)

loss = model.loss
plt.plot(np.arange(len(loss)), np.log10(loss))

plt.show()

A: 932.1335010668406
B: 152.13348416289668

数学的に解いた値と一緒になりました。損失が減っていく様子もみることができますね。

確率的勾配降下法

最急降下法は、サンプル全てを使って勾配を計算しました。確率的勾配降下法では、シャッフルしたデータを1つ〜複数使ってパラメータを更新します。複数使う方法をミニバッチ学習とも言います。

こうすると何がいいかという話なんですが、損失関数は形によって極小値を複数持つことがあります。地形を思い出して欲しいんですが、見えてる谷が一番低いとは限らず、山を越えたらもっと深い谷が見つかることもあります。この一見正しそうだけどそうではない谷底を局所解、もっと低い谷底のことを大域解と呼びます。

最急降下法は初期値の選び方で局所解におちいり易いのですが、確率的勾配降下法はパラメータ更新のたびに勾配計算に使うデータを選び直すので、低い山を越えて大域解にたどり着く可能性が高くなるそうです。

さらに、使うデータが決まっている必要がないのでリアルタイムの処理に向いてるとのことです。

pythonでの実装

先ほどと同じようにStochasticGradientDescentクラスを作りました。中身はパラメータの更新部分以外はほとんど最急降下法と同じです。ミニバッチのサイズは全データに対してどれくらい使うか指定できるようにしました。損失が最低だったパラメータを記憶しておき、それを返すようにしています。

class StochasticGradientDescent:
  def __init__(self, eta=0.01, n_iter=2000, sample_rate=0.1):
    self.eta = eta
    self.n_iter = n_iter
    self.sample_rate = sample_rate
    self.grad = np.zeros((2,))
    self.loss = np.array([])

  def fit(self, X, Y, w0):
    self.w = w0
    self.min_w = w0
    n_samples = int(np.ceil(len(X)*self.sample_rate))
    min_loss = 10**18
    for _ in range(self.n_iter):
      loss = np.sum((Y-self.w[0]-self.w[1]*X)**2)
      if min_loss>loss:
        min_loss = loss
        self.min_w = self.w
      self.loss = np.append(self.loss, loss)
      for i in range(n_samples):
        index = np.random.randint(0, len(X))
        batch_x = X[index]
        batch_y = Y[index]
        self.grad[0] = np.sum(-2*(batch_y-self.w[0]-self.w[1]*batch_x))
        self.grad[1] = np.sum(-2*batch_x*(batch_y-self.w[0]-self.w[1]*batch_x))
        self.w -= self.eta * self.grad

  def predict(self, x):
    return (self.w[0] + self.w[1]*x)

  @property
  def coef_(self):
    return self.min_w[1]

  @property
  def intercept_(self):
    return self.min_w[0]

  @property
  def loss_(self):
    return self.loss

これを実行して、損失の変化も合わせて見てみる。

w0 = np.array([1.,1.])

model = StochasticGradientDescent()
model.fit(x, y, w0)

print("A: ", model.coef_)
print("B: ", model.intercept_)

loss = model.loss
plt.plot(np.arange(len(loss)), np.log10(loss))

plt.show()

A:  925.5203159922469
B:  146.8724188770836

少し厳密解とは異なるのと、少し損失が暴れる感じですね。学習率を低くするとスパイクの高さが低くなるのと、ランダムサンプリングしている関係か、もとのサンプルの相関が低いと厳密解に収束しないのかなと思いました。この辺のパラメータのチューニングは課題ですね。

まとめ

勾配降下法として、最急降下法と買う率的勾配降下法について調べてみました。scikit-learnやディープラーニングの世界ではこれ以外でも様々な勾配方があるんですが、基本は一緒です。こうして順を追って見ていくことで、パラメータのチューニングについても想像がつきやすくなるかもしれません。

次回は回帰の総まとめと過学習や正則化あたりをまとめてみたいと思います。