ニュートン法を Python で解く

# はじめに

ニュートン法は、
$$f(x)=0$$
を $x$ について解くための反復アルゴリズムです。あるいは、
$$f'(x)=0$$
を解くと極値を求めることができます。今回は関数の極大値を求めてみましょう。

# 1 変数関数を解く

ニュートン法のアルゴリズムは、
$$x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}$$
です。今回は、
$$f(x)=\frac{\log x}{x^2}$$
の $1<x<5$ における極大値を求めてみましょう。初期値は $x_0=2$ です。
$$f'(x)=\frac{1-2\log x}{x^3}$$
$$f''(x)=\frac{-5+6\log x}{x^4}$$

## グラフを見る

```py
import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return np.log(x) / x**2

x = np.linspace(1, 5, 100)
y = f(x)

plt.grid(True)
plt.plot(x, y)
plt.show()
```

![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2986389/520d64b2-7765-0203-1e50-d237a268406e.png)



## NumPy を使う

1, 2 階の導関数を自力で求めて計算します。

```py
import numpy as np

def df(x):
    return (1 - 2*np.log(x)) / x**3

def ddf(x):
    return (-5 + 6*np.log(x)) / x**4

# 初期値
x = 2

# 要求精度
epsilon = 1e-5

print("n\tx\tdf(x)")
n = 1

while abs(df(x)) > epsilon:
    # xの更新
    x = x - df(x)/ddf(x)

    print("{}\t{:.5f}\t{:.5f}".format(n, x, df(x)))
    n += 1

print("x = %f" % x)
```

結果

```
n       x       df(x)
1       1.08147 0.66675
2       1.28281 0.23775
3       1.46646 0.07429
4       1.59358 0.01681
5       1.64277 0.00163
6       1.64865 0.00002
7       1.64872 0.00000
x = 1.648721
```

## PyTorch を使う

PyTorch を使えば、微分値の計算を自動でやってくれます。

```py
import torch
from torch.func import hessian, jacfwd

def f(x):
    return torch.log(x) / x**2

# 初期値
x = 2
x = torch.tensor(x, dtype=torch.float64)

# 要求誤差
epsilon = 1e-5

print("n\tx\tdf(x)")
n = 1

df = jacfwd(f)
ddf = hessian(f)

while abs(df(x)) > epsilon:
    # xの更新
    x = x - df(x)/ddf(x)
    
    print("{}\t{:.5f}\t{:.5f}".format(n, x.item(), df(x).item()))
    n += 1

print("x = %f" % x.item())
```

結果

```
n       x       df(x)
1       1.08147 -0.04829
2       1.28281 0.66675
3       1.46646 0.23775
4       1.59358 0.07429
5       1.64277 0.01681
6       1.64865 0.00163
7       1.64872 0.00002
8       1.64872 0.00000
x = 1.648721
```

関数名 | 説明
---|---
`torch.tensor()` | テンソルの作成

## SciPy を使う

SciPy を使えばもっと簡単に解けます。（本当は `scipy.optimize.minimize()` を使うべきですが。）

```py
import numpy as np
from scipy.optimize import newton, approx_fprime

def f(x):
    return np.log(x) / x**2

def df(x):
    return approx_fprime(x, f)

x = newton(df, 2)
print("x = %f" % x.item())
```

結果

```
x = 1.648721
```

# 多変数関数を解く

アルゴリズムは、
$$x_{n+1}=x_n-H^{-1}J$$
です。
今回は、
$$ f(x_0,x_1) = \frac{\log x_0}{x_0^2} - \frac{x_1^2}{37} + \frac{x_1}{7}$$
の $x_0=2, x_1=2$ 付近の極値を求めてみましょう。

## グラフを見る

```py
import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return np.log(x[0])/x[0]**2 - x[1]**2 / 37 + x[1]/7

x0 = np.arange(1, 5, 0.1)
x1 = np.arange(1, 5, 0.1)
x0, x1 = np.meshgrid(x0, x1)

y = f([x0, x1])

plt.contourf(x0, x1, y, cmap='Blues')
plt.show()
```

![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2986389/55172cad-722c-f5be-5434-4eaf67fa5c7e.png)


## PyTorch を使う

```py
import torch
from torch.func import hessian, jacfwd
from torch.linalg import solve

def f(x):
    return torch.log(x[0])/x[0]**2 - x[1]**2/37 + x[1]/7

# 初期値
x = [2, 2]
x = torch.tensor(x, dtype=torch.float64)

# 要求誤差
epsilon = 1e-5

print("n\tx0\tx1\tdf0(x)\tdf1(x)")
n = 1

J = jacfwd(f)
H = hessian(f)

while J(x).sum().abs() > epsilon:
    # xの更新
    x = x - solve(H(x), J(x))

    print("{}\t{:.5f}\t{:.5f}\t{:.5f}\t{:.5f}".format(n,
        x[0].item(), x[1].item(), J(x)[0].item(), J(x)[1].item()))
    n += 1

print("x0 = %f, x1 = %f" % (x[0].item(), x[1].item()))
```

結果

```
n       x0      x1      df0(x)  df1(x)
1       1.08147 2.64286 0.66675 0.00000
2       1.28281 2.64286 0.23775 0.00000
3       1.46646 2.64286 0.07429 0.00000
4       1.59358 2.64286 0.01681 0.00000
5       1.64277 2.64286 0.00163 0.00000
6       1.64865 2.64286 0.00002 0.00000
7       1.64872 2.64286 0.00000 0.00000
x0 = 1.648721, x1 = 2.642857
```

## SciPy を使う

```py
import numpy as np
from scipy.optimize import newton, approx_fprime

def f(x):
    return np.log(x[0])/x[0]**2 - (x[1]**2)/37 + x[1]/7

def df(x):
    return approx_fprime(x, f)

x = newton(df, [2, 2])
print("x0 = %f, x1 = %f" % (x[0].item(), x[1].item()))
```

結果

```
x0 = 1.648721, x1 = 2.642857
```

# 参考

https://zenn.dev/kaitombow/articles/1ee53acb89eee8

https://qiita.com/nariaki3551/items/1ad826df5cf7d0e4ea29

https://qiita.com/tmasada/items/9dee38e5bc1482217493

https://www.hellocybernetics.tech/entry/2020/11/01/155410