ニュートン法の実装と可視化（Python）

#はじめに
たまたま，ニュートン法を後輩に教えるために作っていた限定公開の記事があったので公開しておきます．
ソースコードは、[GitHub](https://github.com/yasu-kun/Newton-s-method)で公開してあります．

#準備＆実行環境
・Mac OSX(windows 10でもOK)
・Python 3.6.7
・matplotlib 3.0.2
・sympy 1.3
・numpy 1.15.4

#実装
１変数の場合のみを、実装していきます。

####ライブラリのインポート


```python:newton_class.py
import sympy as sym
from sympy import sympify
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
```
一通り、今回使用するライブラリです。
[sympy](https://note.nkmk.me/python-sympy-factorization-solve-equation/)は、代数計算のためのライブラリです。
微分や積分、因数分解などいろいろしてくれる便利なライブラリです！
自動微分のプログラムを作ろうと、いろいろ試みたのですが沼にハマりそうだったので、脱却しました。

---
####クラスの定義
```python:newton_class.py
class Newton(object):
    def __init__(self, function):
        self.x = sym.Symbol('x')
        self.function = sympify(function)
        self.y_plot = sym.lambdify(self.x, sympify(function), "numpy")
```
手始めに、コンストラクタを記述します。
ちょっと気持ち悪いことしてます。
・```self.x = sym.Symbol('x')```
ここでは、xは変数だよとシンボル化しています。
・```self.y_plot = sym.lambdify(self.x, sympify(function), "numpy")```
後に関数にndarrayを入れていくので、numpy用の関数も定義してあげました。

---
####Newton法の実装
```python:newton_class.py
def solve(self,seed,iteration=10,verbose=0):
        function = self.function
        x = self.x
        y_plot = self.y_plot

        def cal_diff(diff_function):
            diff_function = sym.diff(diff_function)
            return diff_function
```
・```diff_function = sym.diff(diff_function)```
ここでは、関数を微分しています。
なぜ微分をするかというと、ニュートン法の肝である

```math
x_{n+1} = x_{n} + \frac{f(x_{n})}{f'(x_{n})}\tag{1}
```
を求めていく際に$f'(x_{n})$が必要だからです。

・```sym.diff(diff_function)```
のように、微分したい関数を入れるだけです。
上記のコードを他に代入し、もう一度実行させると２階微分ができます。
（defの中にdefを入れて微分させています。インテリではないです。もう少し良い方法ないでしょうか。）

*iterationとverboseだけ、デフォルトを決めています。seedは勉強のために自分で入れるようにしています。


```python:newton_class.py
#Gitには、ここにplotのコードも入ってます。見やすいように省略しています。
ims = []
for n_iter in range(1,iteration):
    if n_iter == 1:
        Xn_1 = float(seed - function.subs(x,seed) / cal_diff(function).subs(x,seed))

        if verbose==1:
            print("-----ITERATION #{}-----".format(n_iter))
            print("X = {}".format(Xn_1))

            print("ERROR = {}".format(abs(seed - Xn_1)))
        else:
            pass
    else:
        Xn = Xn_1
        Xn_1 = float(Xn - function.subs(x,Xn) / cal_diff(function).subs(x,Xn))
        
        if verbose==1:
            print("-----ITERATION #{}-----".format(n_iter))
            print("X = {}".format(Xn_1))

            print("ERROR = {}".format(abs(Xn - Xn_1)))
        else:
            pass

```
まず、ここでの繰り返しforは、式(1)を繰り返していくための繰り返しです。
そして、一つ目のif文は１回目の計算かそれ以降の計算かを分けています。
最初は初期の計算値が必要だからです。（これも、うまくif文を使わない方法はないでしょうか。）

verboseのif文は、単にロマンがあり導入しました。

```python:newton_class.py
def solve(self,seed,iteration=10,verbose=0):
    #上記のforが入ります。
    return Xn_1, abs(Xn - Xn_1)
```
最後に、近似的な$x$とその時の誤差を出力するようにしました。

---
####実装したコードを実行
```python:newton_class.py
def run():
    y = input('input: ')
    newton = Newton(function=y)
    x,e = newton.solve(iteration=7,seed=10,verbose=1)
    print(x,e)

if __name__ == "__main__":
    run()
```
importすることを考えて、ちょっとお堅く書きました。

・```y = input('input: ')```
関数を入力します。（変数は'x'です（シンボル））
・```x,e = newton.solve(iteration=7,seed=10,verbose=1)```

ここで、いろいろパラメータを設定します。以下パラメータの説明です。
・```iteration```: （近似のための）計算の回数（default=10）
・```seed```　　　: 初期値の設定です。（defaultなし）
・```verbose```  　: 計算過程を表示するかどうか(0:表示しない、1:表示する)、（default=0)

##実行例
$f(x)=x^2 - 5$
上記の式の$f(x)=0$である$x$を求めます。
なぜ、この関数を例として使うかというと、、、
$0=x^2 - 5$
$→x^2 = 5$
$→x = \sqrt{5}$
つまり、$\sqrt{5}$を近似的に求めることができます！

```text:output
input: x**2 - 5
-----ITERATION #1-----
X = 5.25
ERROR = 4.75
-----ITERATION #2-----
X = 3.1011904761904763
ERROR = 2.1488095238095237
-----ITERATION #3-----
X = 2.3567372726441826
ERROR = 0.7444532035462936
-----ITERATION #4-----
X = 2.239157222737191
ERROR = 0.11758004990699167
-----ITERATION #5-----
X = 2.2360701085328496
ERROR = 0.003087114204341379
-----ITERATION #6-----
X = 2.2360679775008054
ERROR = 2.1310320441614294e-06
2.2360679775008054 2.1310320441614294e-06
```
一番最後に```2.2360679```とあります。うまく近似できてます。
###matplotlibによる近似の可視化
どのように、$\sqrt{5}  (=2.23606...)$に近づいてるか、gifで見てみます。
(このplot機能はGitHubにあるコードには実装されています。)
![anim_gif_roop_inf.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/262828/b277449c-7f0a-e1e6-49d1-172b5ea5cd2b.gif)

#新規性
・微分を自動化しました。
・近似の様子をplotしてわかりやすくしました。

##今後の課題

まず、収束しない可能性があります。
つまり、初期値によって結果が変わるので、初期値を何度か試すことが必要である。
もちろん、$f(x)=0$の解が複数する可能性があるのでこちらも、関数が未知の場合は初期値を変えて$x$を探索する必要があります。

また微分不可能な状態の関数に対してのエラー処理はしておりません。


#まとめ
代数計算を汎用性を保ちながら、コーディングする難しさを知りました。
アルゴリズムは1度，1から実装しながら勉強することをお勧めします．

#参考資料
[機械学習のためのニュートン法 (1 変数から多変数まで) - Qiita](https://qiita.com/utahkaA/items/45ebadd70c3020e8c320)