Help us understand the problem. What is going on with this article?

二分法・ニュートン法・勾配降下法をPythonで理解する

More than 1 year has passed since last update.

方程式 f(x) = 0 の解を求める方法として代表的なものに、二分法とニュートン法があります。それをPythonで実装しながら理解していきましょう。

例として、こんな関数を使います。

import math
def f(x):
    return math.exp(x) - 3 * x

ちなみに、以下のように書いても同じ意味になります。

import math
f = lambda x: math.exp(x) - 3 * x

Scipy を使う方法

scipy という科学技術計算用のライブラリが便利です。

from scipy import optimize

二分法

二分法は、関数 y=f(x)が区間[a,b]で連続で、解が1つ存在するときに使えます。たとえば区間[0,1]の中に解があるとわかってる場合は

optimize.bisect(f, 0, 1)
0.619061286737633

あら簡単。

ニュートン法

ニュートン法は、関数y=f(x)が単調連続で変曲点がなく、かつf(x)の導関数が求められるときに使えます。f(x)=0のときのxを求めたいときは

optimize.newton(f, 0)
0.6190612867359452

あら簡単。

もっと詳しい使い方を知りたい

「?」を使うと、パラメーターの説明や使用例などが確認できます。英語ですけど。

optimize.bisect?
optimize.newton?

他にどんな最適化手法があるか知りたい

「optimize?」とすると、.bisect や .newton だけでなく、他の方法についてもいろいろ説明してあります。英語ですけど。

optimize?

ググってももちろん良いんですが、「?」を使うとダイレクトにその説明に行けて、しかもオフラインでも使えるのが良いですね。

Scipy を使わずに実装してみよう

二分法

二分法の理屈は、以下の通りです。

  • 区間[a,b]を2つに分け、その中点をcとする。
  • f(a)*f(c)が負であれば、区間[a,c]の中にx軸との交点があるはず。
  • f(a)*f(c)が正であれば、区間[c,b]の中にx軸との交点があるはず。
  • というのを繰り返す。
import math

def f(x):
    return math.exp(x) - 3 * x

def fd(x):
    return math.exp(x) - 3

eps1 = 0.00001
eps2 = 0.00001

a = 0
b = 1

print("n\tc\tf(c)")
n = 1
while True:
    c = (a + b)/2
    print("{}\t{:.5f}\t{:.5f}".format(n, c, f(c)))
    if f(a)*f(c) < 0:
        b = c
    else:
        a = c
    n += 1
    if abs(f(c)) < eps1 or abs(a - b) < eps2:
        break

print("x = %f" % c)
n   c   f(c)
1   0.50000 0.14872
2   0.75000 -0.13300
3   0.62500 -0.00675
4   0.56250 0.06755
5   0.59375 0.02952
6   0.60938 0.01116
7   0.61719 0.00214
8   0.62109 -0.00232
9   0.61914 -0.00009
10  0.61816 0.00103
11  0.61865 0.00047
12  0.61890 0.00019
13  0.61902 0.00005
14  0.61908 -0.00002
15  0.61905 0.00001
16  0.61906 -0.00000
x = 0.619064

ニュートン法

ニュートン法の理屈は次の通りです。

  • あるxにおける接線を求め、そのx軸との交点を求める。
  • そのx軸との交点のx座標を新たなxとし、そのxにおける接線を求める。
  • というのを繰り返す。
import math

def f(x):
    return math.exp(x) - 3 * x

def df(x):
    return math.exp(x) - 3

eps1 = 1e-5

x = 0 # 初期値

print("n\tx\tf(x)\tdf(x)\tdx")
n = 0
while True:
    print("{}\t{:.4f}\t{:.4f}\t{:.4f}\t{:.4f}".format(n, x, f(x), df(x), f(x)/df(x)))
    x -= f(x)/df(x)
    if abs(f(x)) < eps1:
        break
    n += 1

print("x = %f" % x)
n   x   f(x)    df(x)   dx
0   0.0000  1.0000  -2.0000 -0.5000
1   0.5000  0.1487  -1.3513 -0.1101
2   0.6101  0.0104  -1.1595 -0.0089
3   0.6190  0.0001  -1.1429 -0.0001
x = 0.619061

ニュートン法を実行するクラスを作ってみました

class Newton:
    def __init__(self, f, df):
        self.f = f
        self.df = df

    def _update(self, x):
        return x - f(x)/df(x)

    def solve(self, init_x, n_iter, tol):
        print("n\tx\tf(x)\tdf(x)\tdx")
        x = init_x
        for i in range(n_iter):
            new_x = self._update(x)
            error = abs(new_x - x)
            print("{}\t{:.4f}\t{:.4f}\t{:.4f}\t{:.4f}".format(i, x, f(x), df(x), f(x)/df(x)))
            x = new_x
            if error < tol:
                break
        return x
f = lambda x: math.exp(x) - 3 * x
df = lambda x: math.exp(x) - 3

newton = Newton(f=f, df=df)
result_x = newton.solve(init_x=0, n_iter=100, tol=1e-5)
print("x = {0:.6f}".format(result_x))
n   x   f(x)    df(x)   dx
0   0.0000  1.0000  -2.0000 -0.5000
1   0.5000  0.1487  -1.3513 -0.1101
2   0.6101  0.0104  -1.1595 -0.0089
3   0.6190  0.0001  -1.1429 -0.0001
4   0.6191  0.0000  -1.1428 -0.0000
x = 0.619061

収束条件判定

  1. | f(Xn) | < eps
  2. | Xn+1 - Xn | < eps
  3. | (Xn+1 - Xn) / Xn | < eps
  4. | f(Xn+1) - f(Xn) | < eps

など。

勾配降下法 (Gradient descent)

ついでに勾配降下法についても説明します。こちらは関数の極小値の近似解を求める方法のひとつですが、計算値と実際の値との誤差をどんどん小さくしていこうという意味では似たようなものです。

数値微分(微分の数値的近似, 1変数の場合)

def numerical_differentiation(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2 * h)
def f(x):
    return x ** 2 + 3 * x + 1
f(0.5)
2.75

偏微分の数値的近似, 複数変数の場合

import numpy as np
def numerical_partial_differential(f, x):
    h = 1e-4
    grad = np.zeros_like(x)
    for idx_i in range(len(x)):
            tmp_val = x[idx_i]

            # f(x+h)
            x[idx_i] = tmp_val + h
            fxh1 = f(x)

            # f(x-h)
            x[idx_i] = tmp_val - h
            fxh2 = f(x)

            grad[idx_i] = (fxh1 - fxh2) / (2 * h)
            x[idx_i] = tmp_val
    return grad
import math
def f(x):
    return (x[0] + math.pi)**2 + (3 * x[1] - math.sqrt(2))**2
x = np.array([0.0, 0.0])
numerical_partial_differential(f, x)
array([ 6.28318531, -8.48528137])

勾配降下法

def gradient_decent(f, init_x, lr=0.1, step_num=100):
    x = init_x
    for i in range(step_num):
        grad = numerical_partial_differential(f, x)
        x -= lr * grad
    return x
init_x = np.array([5.0, 5.0])
gradient_decent(f, init_x)
array([-3.14159265,  0.47140452])

どのように最適化が進んだか履歴を追ってみましょう。

import copy
def gradient_decent(f, init_x, lr=0.1, step_num=100):
    x = init_x
    history = []
    for i in range(step_num):
        history.append(copy.deepcopy(x))
        grad = numerical_partial_differential(f, x)
        x -= lr * grad
    return x, np.array(history)
init_x = np.array([5.0, 5.0])
result_x, history = gradient_decent(f, init_x)
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

plt.plot(history[:, 0], history[:, 1])
plt.grid()
plt.show()

output_14_0.png

課題

以下の問題を、scipy を使わずに解いてください(scipyを使って解いた結果との比較はOKです)。また、pi は円周率とします。

  • f(x) = x2 - 2 x + 1 のとき、y = f(x) のグラフを描いてください。また f(x) = 0 の近似解をニュートン法により求めてください。そのときに、様々な eps や、収束条件1〜4を試してみて、収束状況の違いを確認してください。

  • f(x) = ex sin x - 2 x (0 < x < pi) のとき、y = f(x) のグラフを描いてください。また f(x) = 0 の近似解をニュートン法により求めてください。

  • f(x) = tanh x + x + 2 のとき、y = f(x) のグラフを描いてください。また f(x) = 0 の近似解を2分法により求めてください。

  • 二分法を実行するクラスを作ってください。またその動作確認をしてください。

  • y = 2 * sin (x1 - pi) + 3 * (x2 + pi)2 に対して勾配降下法を適用し、極小値を求めてください。また、学習率 lr や初期値を変えてその挙動を確認してください。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away