最急降下法とニュートン法の比較

autogradを使って、最急降下法（steepest descent method）とニュートン法（Newton's method）の比較をしてみた。

最急降下法のおさらい

1. 適当な初期点xを選ぶ
2. xでの勾配ベクトルを求めて勾配ベクトル方向（最も急な方向）に点を移動させる
3. 収束するまで2.を繰り返す

ニュートン法のおさらい
関数fの点x_nのまわりでのテイラー展開は以下のようになる。

ただし、3次以上の項は無視している。

最適化問題で求めたいのはfの停留点である。

ということで、上のテイラー展開を使って微分が零になる点を近似的に計算してみる。

なので、fをxで微分したものと、fをΔxで微分したものは同じ。

よって、

となる。この値が0になるようなΔxを求めると、

である。

つまりこの方向に進めば、fの微分が0になりそうな場所に向かって、xを逐次的に更新することができる。

ニュートン法のアルゴリズムを以下に示す。

1. 適当な初期点xを選ぶ
2. 上記の手順でx = x + Δxと更新する
3. 収束するまで2.を繰り返す

実験
autogradを使って、最急降下法とニュートン法を実装する。
必要なライブラリのimport。

import autograd
import matplotlib.pyplot as plt
import numpy as np

最急降下法の実装。
グラフに描画したいので、xの最適値だけじゃなくて、xの更新履歴も返すようにしている。

def steepest_descent(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x -= alpha * g(x)
    history.append(np.array(x))
  return x, history

ニュートン法の実装。

def newton(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    H = autograd.hessian(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x += alpha * np.linalg.solve(H(x), -g(x))
    history.append(np.array(x))
  return x, history

軸がずれた楕円形の関数を使って、両手法を比較してみる。

def f(x, y):
  s = 2*x + y
  t =  -x + y
  return 4 * s**2 + 3 * t**2

# draw contour of f(x, y)
delta = 0.01
x = np.arange(-2.2, 0.2, delta)
y = np.arange(-2.2, 0.2, delta)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
plt.contour(X, Y, Z, 80, colors='black')

# execute steepest descent
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = steepest_descent(wrapper, x0, alpha=0.001)
for x in history:
  plt.plot(x[0], x[1], 'ro', markeredgewidth=0.0)

# execute newton method
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = newton(wrapper, x0)
for x in history:
  plt.plot(x[0], x[1], 'bo', markeredgewidth=0.0)
  
plt.show()

結果は以下のとおり。
最急降下法（赤点）の場合は、等高線と直交する方向にxが更新されているのが分かる。
それに対して、ニュートン法（青点）の場合は、近似的に計算した関数の停留点に向かって真直ぐxが更新されているのが分かる。