方程式 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

収束条件判定

| f(X_n) | < eps
| X_n+1 - X_n | < eps
| (X_n+1 - X_n) / X_n | < eps
| f(X_n+1) - f(X_n) | < 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()

課題

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

f(x) = x² - 2 x + 1 のとき、y = f(x) のグラフを描いてください。また f(x) = 0 の近似解をニュートン法により求めてください。そのときに、様々な eps や、収束条件1〜4を試してみて、収束状況の違いを確認してください。
f(x) = e^x 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 の近似解を２分法により求めてください。
二分法を実行するクラスを作ってください。またその動作確認をしてください。
y = 2 * sin (x₁ - pi) + 3 * (x₂ + pi)² に対して勾配降下法を適用し、極小値を求めてください。また、学習率 lr や初期値を変えてその挙動を確認してください。