TOC

7. Newton 法と勾配降下法

7.1 Newton 法アルゴリズム

\(k=0\) とおく,
初期値 \(x_0\) を設定する.
以下を \(x\) の変化の絶対値 \(\epsilon = \lvert x_{k+1} - x_k \rvert\) が許容範囲 \(\epsilon_0\) 以下になるまで繰り返す.
- \(x_{k+1} = x_k - f(x_k) / f^\prime(x_k)\) を算出して \(x\) を更新する.
- \(x\) の変化の絶対値 \(\epsilon = \lvert x_{k+1} - x_k \rvert\) を計算する.
- \(k\) を 1 増やす.

注意

極値を求めるには \(f^\prime(x) = 0\) を解く.
極値（極大値か極小値か, どちらでもないか）の判別が必要.

7.2 勾配降下法（最急降下法, gradient descent）アルゴリズム

\(f^\prime(x)\) を算出する.
次のように \(x\) の推定値を更新する.

\[x_{k+1} = x_k - \gamma f^\prime(x_k)\]

注意

\(f^\prime(x)\) が十分小さくなったらそのときの \(x\) を極小値とする.
極大値は \(y = -f(x)\) の極小値を求める.
極値かどうかの判定が必要
\(\gamma\) はステップ幅（step size）.

例題 1

次の関数の極大値と極小値および極値を与える \(x\) の値を Newton 法と勾配降下法により、それぞれ求めよ.

\[y = x^4 - x^2\]
反復回数 5 回までの値を求めよ. 収束判定のプログラミングは行わなくてよい.
勾配降下法では, ステップ幅の影響

[1]:

import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline


def func(x):
    return x ** 4 - x ** 2


def dfunc(x):
    return 4 * x ** 3 - 2 * x


def ddfunc(x):
    return 12 * x ** 2 - 2


fig, ax = plt.subplots()
xmin, xmax = -1, 1
xs = np.linspace(xmin, xmax)
ys = func(xs)
dy = dfunc(xs)
ax.plot(xs, ys)
ax.plot(xs, dy)
# ax.plot(xs, ddfunc(xs))
ax.plot([xmin, xmax], [0, 0], lw=.5, c='r')
ax.set_ylim(-2, 2)

[1]:

(-2.0, 2.0)

../_images/exercise_process_engineering_2021_exercises_070_4_1.png

[2]:

# Newton's method

x_0 = 1 # 初期値
x = x_0
x_p = x - dfunc(x) / ddfunc(x)
print(x, x_p)
x = x_p
x_p = x - dfunc(x) / ddfunc(x)
print(x, x_p)
x = x_p
x_p = x - dfunc(x) / ddfunc(x)
print(x, x_p)
x = x_p
x_p = x - dfunc(x) / ddfunc(x)
print(x, x_p)
x = x_p
x_p = x - dfunc(x) / ddfunc(x)
print(x, x_p)

1 0.8
0.8 0.7211267605633803
0.7211267605633803 0.7075053183719765
0.7075053183719765 0.7071071176773481
0.7071071176773481 0.7071067811867877

[22]:

# Gradient descent

gmm = .25
y = -1
y_p = y - gmm * dfunc(y)
print(y, y_p)
y = y_p
y_p = y - gmm * dfunc(y)
y = y_p
print(y, y_p)
y = y_p
y_p = y - gmm * dfunc(y)
print(y, y_p)
y = y_p
y_p = y - gmm * dfunc(y)
print(y, y_p)
y = y_p
y_p = y - gmm * dfunc(y)
print(y, y_p)

-1 -0.5
-0.625 -0.625
-0.625 -0.693359375
-0.693359375 -0.7067084684967995
-0.7067084684967995 -0.707106444695907

例題 2

上記の問題を, scipy.optimize.newton を用いて解け.

[25]:

import scipy.optimize

x_0 = .5

res = scipy.optimize.newton(dfunc,-1, fprime=ddfunc)
print(res)

-0.7071067811865476

例題 3

次の関数について以下の問に答えよ.

\[f(x, y) = x^2 + y^2 - 2y\]

1. 勾配とヘッセ行列を求めよ.
1. 停留点を求め, 極小値あるいは極大値を与えるかどうかを判定せよ.
1. 極大値か極小値がある場合, 数学ツールでそれを求めよ.
1. 2 次元の Newton 法で, 極小値を与える \((x, y)\) を求めよ

【参考】多変数関数 \(f:\mathbb R^n \to \mathbb R\) の \(\operatorname{grad} f (\boldsymbol x) = 0\) の解（停留点）を求める場合

グラディエント（勾配）とヘッセ行列の定義

\begin{equation} \nabla f (\boldsymbol x) = \operatorname{grad} f(\boldsymbol x) = \begin{bmatrix}\frac{\partial f}{\partial x_1} \\ \vdots \\ \frac{\partial f}{\partial x_n} \end{bmatrix} \tag{3.7} \end{equation}

\begin{equation} H f(\boldsymbol x) = \begin{bmatrix} \frac{\partial^2 f}{\partial x^2} & \dots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \dots & \frac{\partial^2 f}{\partial x_n^2} \end{bmatrix} \tag{3.8} \end{equation}

\(\operatorname{grad} f (\boldsymbol x) = 0\) の解を求めるアルゴリズム

\[\boldsymbol x^{[k+1]} = \boldsymbol x^{[k]} - \underbrace{\left[ H f(\boldsymbol x^{[k]}) \right]^{-1}}_{\mathbb R^n \times \mathbb R^n} \underbrace{ \operatorname{grad} f(\boldsymbol x^{[k]})}_{\mathbb R^n} \tag{3.9}\]

解答例 3

\[\begin{split}\operatorname{grad} f = \begin{bmatrix} 2x \\ 2 y -2 \end{bmatrix}\end{split}\]

\[\begin{split}Hf = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}\end{split}\]

\[2x = 0, \quad 2 y - 2 = 0\]

停留点は上の連立方程式の解だから, \((x, y) = (0, 1)\)
- この点におけるヘッセ行列の固有値は \(\lambda = 2, 2\),
- これは重複しているが縮退はしていない（固有ベクトルは 2 つある.）
- よって, 停留点は極小値を与える.

[77]:

#### 例題 3 (3)
import numpy as np


def fn(x):
    return x[0]**2 + x[1]**2 - 2 * x[1]

h = np.array([[2, 0], [0, 2]])
eigval, eigvec = np.linalg.eig(h)
print(eigval)
print(eigvec)

x_0 = [1, 1]
res = scipy.optimize.minimize(fn, x_0)
print(res)

[2. 2.]
[[1. 0.]
 [0. 1.]]
      fun: -0.9999999999999997
 hess_inv: array([[5.00000005e-01, 7.37681254e-09],
       [7.37681254e-09, 1.00000000e+00]])
      jac: array([-7.45058060e-09, -2.23517416e-08])
  message: 'Optimization terminated successfully.'
     nfev: 9
      nit: 2
     njev: 3
   status: 0
  success: True
        x: array([-9.44232017e-09,  9.99999985e-01])

[35]:

#### 例題 3 (4) (Newton in 2d)
import numpy as np


def fn(x):
    return x[0]**2 + x[1]**2 - 2 * x[1]


def grad_fn(x):
    return np.array([2 * x[0], 2 * x[1] - 2])


def hess_fn(x):
    return np.array([[2, 0], [0, 2]])


x_0 = np.array([2, 5])
x = x_0
count = 0
print(x)
while count < 5:
    count +=1
    hinv = np.linalg.inv(hess_fn(x))
    grad = grad_fn(x)
    dx = np.matmul(hinv, grad)
    x = x - dx
    print(x)

[2 5]
[0. 1.]
[0. 1.]
[0. 1.]
[0. 1.]
[0. 1.]

演習問題

問題 1

数学ツールを使うこと.
適宜グラフを描くこと.
次の関数の極値を与える \(x\) とその極値および極小値か極大値か, を, すべての停留点について調べよ.

\[f(x) = x^5 - 4 x^3 + 2 x\]

問題 2

\[f (x, y) = x^3 + y^3 + x^2y^2 - x - y\]

1. 勾配ベクトルとヘッセ行列が次のようになることを確かめよ.

\[\begin{split}\operatorname{grad} f = \begin{bmatrix} 3x^2 + 2xy^2 - 1 \\ 3 y ^ 2 + 2x^2 y - 1 \end{bmatrix}\end{split}\]

\[\begin{split}H f = \begin{bmatrix}6 x + 2y^2 & 4xy \\ 4xy & 6y + 2x^2 \end{bmatrix}\end{split}\]

1. 関数 \(f\) に対して初期値を適切に設定して Newton 法を適用して, 極小値を求めよ（今回はちゃんと極小値があります）.

◀ Previous

TOC

プロセス工学演習（後半）