ニュートン法（2変数）

ニュートン法（１変数）

ニュートン法（2変数）

2変数のニュートン法も作りました。

凸関数 $f(x, y)$ の最小解と最小値をニュートン法で求める

\begin{aligned}
&\nabla f(x, y)=\left[\begin{array}{ll}
f_{x}(x, y) \\
f_{y}(x, y)
\end{array}\right] \\
&\nabla^{2} f(x, y)=\left[\begin{array}{ll}
f_{x x}(x . y) & f_{x y}(x, y) \\
f_{y x}(x, y) & f_{y y}(x, y)
\end{array}\right] \\
&\left(\begin{array}{l}
x_{new} \\
y_{new}
\end{array}\right)=\left(\begin{array}{l}
x_{old} \\
y_{old}
\end{array}\right)-\nabla^{2} f(x_{old}, y_{old})^{-1} \nabla f(x_{old}, y_{old})
\end{aligned}

更新式から最小解 $(x,y)$ を求めて、最後に最小値を出力する

実装

import matplotlib.pyplot as plt
import numpy as np


# 実装
def newton_method_2val(func0,func1, func2, initial_x, initial_y, EPS):
    """
    func0: f(x,y)
    func1: f_x(x,y)
    func2: f_y(x,y)
    funcの部分はlambda関数で設定する。
        func = lambda x,y : x+y+3
        func(1,3) = 7
    EPS：更新幅を止める基準
    """
    count = 0
    # xの更新値の保存
    x_new_list = [ ]
    y_new_list = [ ]
    diff = 10 # EPS基準で止めに行くが、止まらない値で設定しておく
    
    print("========== processing ==========")
    
    while diff >= EPS:
        count += 1
        f_x = func1(initial_x, initial_y)
        f_y = func2(initial_x, initial_y)

        # nabla^2 f(x,y)
        h = 1e-10 # 微小区間
        # 中心差分公式で微分
        f_xx = (func1(initial_x + h, initial_y) - func1(initial_x - h, initial_y))/(2*h)
        f_xy = (func1(initial_x, initial_y+h) - func1(initial_x , initial_y-h))/(2*h)
        f_yx = (func2(initial_x + h, initial_y) - func2(initial_x - h, initial_y))/(2*h)
        f_yy = (func2(initial_x, initial_y+h) - func2(initial_x , initial_y-h))/(2*h)
        
        # 更新に必要なものを行列形式にまとめる
        initial = np.array([[initial_x], [initial_y]]) # 初期値格納
        nabla = np.array([[f_x], [f_y]]) # nabla(x,y)
        nabla2 = np.array([[f_xx, f_xy],[f_yx, f_yy]]) # nabla^2(x,y)
        inv_nabla2 = np.linalg.inv(nabla2) # nabla^2(x,y) の逆行列
        
        # 更新値を求める
        new = initial - np.dot(inv_nabla2, nabla)
        # count と 更新値
        print("{}回目:x_new = {},y_new={}".format(count,new[0,0],new[1,0]))
        x_new_list.append(new[0,0]) # xの更新値をリストに追加
        y_new_list.append(new[1,0]) # yの更新値をリストに追加
        
        
        # 更新幅
        diff_x = abs(initial_x - new[0,0])
        diff_y = abs(initial_y - new[1,0])
        # diffの大きい方でwhileを止める基準にしたい
        diff = max(diff_x,diff_y)
        # 値の更新
        initial_x = new[0,0]
        initial_y = new[1,0]

        
        if count == 100: # 100回で収束しなければ
            diff = 0 # diffを強制的に0にして
            print("========== Do not CONVERSION =========")  # while 文から抜ける
    
    print("========== Result ==========")
    print("solution:x = {}".format(new[0,0]))
    print("solution:y = {}".format(new[1,0]))
    value = func0(new[0,0],new[1,0]) # 解を代入した値
    print("min_value:{}".format(value))
    
    # グラフ描画
    trial = np.arange(1, count+1, 1) # 1 ~ n回

    fig = plt.figure(figsize=(15, 2.5))
    ax = fig.add_subplot(1, 2, 1)
    x_value = np.array(x_new_list)
    plt.plot(trial, x_value, color="red")
    plt.title("x_process")
    plt.xlabel('trial')
    plt.ylabel('x_new')

    ax = fig.add_subplot(1, 2, 2)
    y_value = np.array(y_new_list)
    plt.plot(trial, y_value, color="blue")
    plt.title("y_process")
    plt.xlabel('trial')
    plt.ylabel('y_new')


    plt.show()

投入

$$f(x, y)=\frac{x^{4}}{2}+\frac{x^{3}}{3}-2 x^{2} y+4 x+3 y^{2}+4 y+4$$

newton_method_2val(lambda x,y: (1/2)*x**4+(1/3)*x**3-2*(x**2)*y+3*y**2+3*x-1*y+4,
                   lambda x,y: 2*x**3+x**2-4*x*y+3, 
                   lambda x,y: 6*y-2*x**2-1,
                   9,3, 0.0001)

========== processing ==========
1回目:x_new = 6.956495810248076,y_new=14.90540384620938
2回目:x_new = 4.61060249609451,y_new=5.418088454657781
3回目:x_new = 3.129220719387754,y_new=2.6991904039382333
4回目:x_new = 2.021337899706909,y_new=1.1194672031808257
5回目:x_new = 1.1684879443600464,y_new=0.3793366434284833
6回目:x_new = 0.3027062364792097,y_new=-0.05264901159486629
7回目:x_new = -2.290661442659664,y_new=-0.3261419277396098
8回目:x_new = -2.414577217675341,y_new=2.1049424354116724
9回目:x_new = -2.5851419887405522,y_new=2.384622254255956
10回目:x_new = -2.5704464991029425,y_new=2.368993071767485
11回目:x_new = -2.5702482687485855,y_new=2.368725374647714
12回目:x_new = -2.5702482387915295,y_new=2.3687253363369853
========== Result ==========
solution:x = -2.5702482387915295
solution:y = 2.3687253363369853
min_value:-4.38238055980867

前の記事では縦軸を差分で出していたので、trialを増すごとに差分が小さくなっていく様子が見れますが、今回は解の更新の推移を表示しているので、必ずしも小さくなるとは限りませんが収束はしていきます。

答え合わせ

答え合わせに使用した wolframAlpha
https://www.wolframalpha.com/