機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

ニュートン法(多変数)の実験

ニュートン法(多変数の場合)で導出した


\mathbf{x} = \mathbf{x}^{\prime} - \bar{\mathbf{H}}^{-1} \nabla \bar{f} \tag{1}

を使って、関数の極値を求めてみたいと思います。

関数は


\displaystyle f(x_1, x_2)=\frac{1}{2 }x_1^{4} - 2x_2^{2} x_2 + 4x_2^{2} + 8x_1 + 8x_2 \tag{2}

としました。*1 この関数の形状と等高線はこちらのようになっています。

f:id:opabinia2:20180807221601p:plain

さて、この関数のヘッセ行列\mathbf{H}と勾配\nabla fは、


\mathbf{H} = 
 \left(
    \begin{array}{cc}
      6x_1^{2}-4x_2  & -4x_1 \\
      -4x_1 & 8 
    \end{array}
  \right) \tag{3}


\nabla f =
 \left(
\begin{array}{c}
2x_1^{3}-4x_1 x_2 + 8 \\
-2x_1^{2}+8x_2 + 8
    \end{array}
  \right)
\tag{3}

です。あとは初期値を決めて式(1)を繰り返すだけです。反復によって解へ近づいていく過程をグラフにしました。 f:id:opabinia2:20180807223301p:plain 初期値は(x_1, x_2)=(3,1)としました。比較として、最急降下法(SD:steepest descent)も載せています。最急降下法は解に向かってジグザグしている感じですが、ニュートン法は無駄な動きがあまりなく収束の速さが見て取れます。真の解との距離(ノルム)をグラフにすると、よりはっきりわかります。 f:id:opabinia2:20180807223842p:plain 横軸が反復回数で、縦軸は真の解\mathbf{\alpha}と近似解の差のノルムです。ニュートン法の速さがはっきりわかりますね。

今回のコードです。

# ニュートン法2変数
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm


# 最適化の対象
def f(x1,x2):
    return 0.5*(x1**4) - 2*(x1**2)*x2 + 4*(x2**2) + 8*x1 +8*x2

# ヘッセ行列
def Hesse(x1,x2):
    return np.array([[6*(x1**2)-4*x2, -4*x1], [-4*x1, 8]])

# 勾配
def nablaf(x1,x2):
    return np.array([2*(x1**3)-4*x1*x2+8, -2*(x1**2) + 8*x2 +8])


x1, x2 = np.meshgrid(np.linspace(-4, 4, 100), np.linspace(-4, 4, 100))
y = np.vectorize(f)(x1, x2)

# 初期値
X = np.array([3.0, 1.0])

point = np.copy(X)
point2 = np.copy(X)
er = []
er2 = []

# 解(http://www.wolframalpha.com/で求めたもの)
A = np.array([-1.36466, -0.534429])

# ニュートン法
for i in range(5):
    H = Hesse(X[0], X[1])
    H_inv = np.linalg.inv(H)
    X -= np.dot(H_inv, nablaf(X[0], X[1]))
    er.append(np.linalg.norm(A-X))
    point = np.vstack([point, X])

# 最急降下法
X = np.array([3.0, 1.0])
for i in range(10):
    X -= 0.1*nablaf(X[0], X[1])
    er2.append(np.linalg.norm(A-X))
    point2 = np.vstack([point2, X])


plt.contour(x1, x2, y, cmap=cm.coolwarm, levels=np.linspace(0,400,80))
plt.plot(point[:, 0], point[:, 1],color="green", marker=".", label="Newton")
plt.plot(point2[:, 0], point2[:, 1],color="black", marker=".", label="SD")
plt.legend(loc="best")
plt.show()


plt.yscale("log")
plt.xlabel("iteration")
plt.plot(np.arange(1, 6), er, color="green",marker=".", label="Newton")
plt.plot(np.arange(1, 11), er2, color="black",marker=".", label="SD")
plt.legend(loc="best")
plt.show()

*1:書籍「基礎系 数学 最適化と変分法」から引用しました