最急降下法を用いた連立方程式の近似解の算出について

はじめに

　連立方程式は、例え線形的なものであっても行列のサイズが大きくなれば手計算で解くのは厳しくなる。そのため、計算機を用いて連立方程式を解くことを考える。一般的に用いられる手法としては、ガウス法（参考文献1）がある。この手法は、行基本変形を用いる最も基本的な方法であり、厳密解を求めたい場合に有効である。ただし、ガウス法は、線形的な方程式のみしか適用できない。
　そこで今回は、連立方程式の近似解を求めるために、最急降下法を用いる。
具体的には、現状の数値解での出力ベクトルと目標ベクトルとの偏差を目的関数とする。
その目的関数を最小にするため、最急降下法を用いる。

問題

以下の連立方程式(参考文献1)の近似解を最急降下法を用いて算出する。

{\begin{align}
\begin{cases}
5x-y-z=0 \\
2x+y-3z = -5 \\
x+y+z=6
\end{cases}
\end{align}
}

この方程式の厳密解は$x=1,y=2,z=3$となる。これは、ガウス法により確かめることができる。

モデル化

上記の連立方程式を一般化し、以下のように表記する。
($A$は3x3の行列であり、$x_{ex},x_{opt},b$はベクトル)

Ax_{ex}=b

ここで、近似値を$x_{opt}$とすると、正の実数で表記される目的関数(誤差)は、以下のように与えられる。

err=|b-Ax_{opt}|

この$err$を最小化する$x_{opt}$を数値解とする。
ただし、最小化のアルゴリズムには、最急降下法（参考文献2）を用いた。

プログラム

問題である連立方程式の数値解を求めるべく、以下のプログラム(Python)を作成した。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#変数をしまう配列
err_ary=[]
x_opt_ary=[]
i_ary=[]
#微小量
h=1e-5

# solve "mat @ x_exact = b" about x_exact
# x_exact=[1,2,3] is the solution
mat = np.array([[5,-1,-1], [2,1,-3], [1,1,1]])
b = np.array([0,-5,6])

#試行回数
num=300
#学習率(最急降下法)
alpha=0.01
#初期値
x_opt=np.random.rand(3)
#目的関数
def calc_x(x_opt):
    err= np.linalg.norm(b-mat@x_opt)
    return err
for i in range(num):
    #目的関数の偏微分を求める。
    x_opt1=x_opt+np.array([h,0,0])
    x_opt2=x_opt+np.array([0,h,0])
    x_opt3=x_opt+np.array([0,0,h])
    
    err1=calc_x(x_opt1)
    err2=calc_x(x_opt2)
    err3=calc_x(x_opt3)
    
    derr_dx=(err1-calc_x(x_opt))/h
    derr_dy=(err2-calc_x(x_opt))/h
    derr_dz=(err3-calc_x(x_opt))/h
    
    #最急降下法の更新式
    x_opt=x_opt-alpha*np.array([derr_dx,derr_dy,derr_dz])
    err_ary.append(calc_x(x_opt))
    x_opt_ary.append(x_opt.copy())
    i_ary.append(i)

#目的関数(誤差)の変化    
plt.plot(i_ary,err_ary)
plt.xlabel("iteration")
plt.ylabel("error") 
plt.savefig("error_連立方程式.png")
plt.show()

#厳密解x_exact=[1,2,3]を用意して、数値解と比較する。    
x_exact = np.array([1,2,3])

#解の比較
colors = ["tab:blue","tab:orange","tab:green"]
for k,name in enumerate(["x","y","z"]):
    plt.plot(i_ary,[x[k] for x in x_opt_ary],color=colors[k],label=f"{name} (数値解)")
    plt.axhline(x_exact[k],color=colors[k],linestyle="--",label=f"{name}={x_exact[k]} (厳密解)")
plt.xlabel("iteration")
plt.ylabel("x_opt VS x_ex")
plt.legend()
plt.savefig("x_opt_連立方程式.png")
plt.show()

結果

上記のプログラムを実行すると、以下のグラフが得られた。

目的関数の推移

このように、試行回数を増加させると、目的関数はほぼ単調に減少し0に漸近した。しかし、約200回以降で目的関数の値が振動してしまった。これは、学習率が大きくこれ以上最適化できなかったためである。

数値解の推移

数値解と厳密解を比較した場合、試行回数が増加するにつれて、厳密解に数値解が漸近した。
ただし、目的関数と同様に、試行回数が約200回以降で振動してしまった。

まとめ

連立方程式は、行列のサイズが大きくなるにつれて解くために必要な計算時間が著しく増加する。
これは、ガウス法を用いた場合、顕著に現れる。また、ガウス法は線形的な方程式にしか適用できない。

そこで、それらの解決策の一つとして、最急降下法を用いた近似解の推定方法を紹介した。
ただし、最急降下法は、学習率を適切に選択しないと、数値解や目的関数が振動してしまう場合がある。また、数値解を求めるためには、適当な初期値が必要になることも忘れてはいけない。この初期値については、別の手法もしくは経験的な手法に頼らざるを得ない。

参考文献

1.連立方程式とガウス法

2.線形回帰と最急降下法