ニュートン法（多変数）の実験 - 機械学習に詳しくなりたいブログ

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

$\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_1^{2} x_2 + 4x_2^{2} + 8x_1 + 8x_2 \tag{2}$

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

さて、この関数のヘッセ行列 $\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)を繰り返すだけです。反復によって解へ近づいていく過程をグラフにしました。初期値は $(x_1, x_2)=(3,1)$ としました。比較として、最急降下法(SD:steepest descent)も載せています。最急降下法は解に向かってジグザグしている感じですが、ニュートン法は無駄な動きがあまりなく収束の速さが見て取れます。真の解との距離（ノルム）をグラフにすると、よりはっきりわかります。横軸が反復回数で、縦軸は真の解 $\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:書籍「基礎系数学最適化と変分法」から引用しました