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

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

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

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

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

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

# Scipy を使う方法

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

```python
from scipy import optimize
```

## 二分法

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

```python
optimize.bisect(f, 0, 1)
```




    0.619061286737633


あら簡単。

## ニュートン法

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

```python
optimize.newton(f, 0)
```

    0.6190612867359452


あら簡単。

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

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

```python
optimize.bisect?
```

```python
optimize.newton?
```

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

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

```python
optimize?
```

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


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

## 二分法

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

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

```python
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における接線を求める。
* というのを繰り返す。


```python
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


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

```python
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
```


```python
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(X<sub>n</sub>) | < eps
2. | X<sub>n+1</sub> - X<sub>n</sub> | < eps
3. | (X<sub>n+1</sub> - X<sub>n</sub>) / X<sub>n</sub> | < eps
4. | f(X<sub>n+1</sub>) - f(X<sub>n</sub>) | < eps

など。


# 勾配降下法 (Gradient descent)

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

## 数値微分（微分の数値的近似, 1変数の場合）


```python
def numerical_differentiation(f, x):
    h = 1e-4
    return (f(x+h) - f(x-h)) / (2 * h)
```


```python
def f(x):
    return x ** 2 + 3 * x + 1
```


```python
f(0.5)
```




    2.75



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


```python
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
```


```python
import math
def f(x):
    return (x[0] + math.pi)**2 + (3 * x[1] - math.sqrt(2))**2
```


```python
x = np.array([0.0, 0.0])
numerical_partial_differential(f, x)
```




    array([ 6.28318531, -8.48528137])



# 勾配降下法


```python
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
```


```python
init_x = np.array([5.0, 5.0])
gradient_decent(f, init_x)
```




    array([-3.14159265,  0.47140452])



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


```python
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)
```


```python
init_x = np.array([5.0, 5.0])
result_x, history = gradient_decent(f, init_x)
```


```python
%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](https://qiita-image-store.s3.amazonaws.com/0/211162/efeea32f-6a9c-82e2-5541-747f4802babf.png)


# 課題

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

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

* f(x) = e<sup>x</sup> 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<sub>1</sub> - pi) + 3 * (x<sub>2</sub> + pi)<sup>2</sup> に対して勾配降下法を適用し、極小値を求めてください。また、学習率 lr や初期値を変えてその挙動を確認してください。