ほぼPython

Not技術ブログBut勉強ブログ 内容には誤りがあることが多いです

Sympyを使って2変数関数の極大・極小・鞍点を求める

お久しぶりです。

 

大学院で最適化問題に関する授業を取っています。完全に専門外の分野ですが、興味があるので勉強しています。

 

最適化問題に入る前に基礎的な数学の復習として2変数関数の極大・極小・鞍点を求めよという課題が出たのでPythonを使って解いてみました。

 

今回はSympyというPythonで記号計算を行うためのライブラリを使用しました。

 

まず、今回計算に使う2変数関数はこちらです。

{}

\[ f(x,y) = 2x^3 + 4xy^2 - 10xy + y^2\]

 

この関数の極大・極小・鞍点を求めていきます。

 

そのためにはまず、停留点(stationary points)を求めなければなりません。

 

停留点を求める

停留点とは f(x,y)xあるいはy偏微分した2つの偏導関数が0となる点のことです。

もう少し一般化すると、実多変数の可微分関数に対してすべての偏微分が0になるような定義域内の値らしいです。

 

つまり、以下の関係を満たす(x,y)の組を求めれば良いです。



\[ f_x(x,y) = f_y(x,y) = 0\]



Sympyを使うとこれを簡単に解くことが出来ます。

from sympy import Symbol, solve, Derivative, Matrix, simplify
x = Symbol('x')
y = Symbol('y')

f = 2*x**3 + 4*x*y**2 - 10*x*y + y**2

f_x = Derivative(f,x).doit() #xで1階偏微分
f_y = Derivative(f,y).doit() #yで1階偏微分

stationary_points = solve([f_x,f_y],[x,y]) #関数の停留点を求める
stationary_points = list(map(lambda x:simplify(x),stationary_points)) #簡単化


停留点を出力すると以下のようになりました。

 [(0, 0), (1, 1), (-3/4 - sqrt(6)/12, -sqrt(6)/8 + 2), (-3/4 + sqrt(6)/12, sqrt(6)/8 + 2)]

さて、これで停留点を求めることが出来ました。

停留点は、極大・極小の候補となります。また、極値ではない停留点は鞍点と呼ばれます。

これらの停留点が、極大点・極小点・鞍点のどれに分類されるか決定するために、ヘッセ行列を使います。


ヘッセ行列とその行列式を求める

ヘッセ行列とは以下の形をした行列のことです。



\[
H(x,y) = \begin{bmatrix}
f_{xx}(x,y) & f_{xy}(x,y) \\
f_{yx}(x,y) & f_{yy}(x,y) \\
\end{bmatrix}
\]



ヘッセ行列の行列式は以下のように表されます。



\[
det H(x,y) = f_{xx}(x,y)*f_{yy}(x,y) - f_{xy}(x,y)*f_{yx}(x,y)
\]



この行列式 det H(x,y)f_{xx}の正負により、基本的には、極大・極小・鞍点の判別が可能となります。

では、ヘッセ行列とヘッセ行列の行列式を求めます。

f_xx = Derivative(f, x, 2).doit() #xで2階偏微分
f_yy = Derivative(f, y, 2).doit() #yで2階偏微分
f_xy = Derivative(f_x,y).doit() #xで偏微分してからyで偏微分
f_yx = Derivative(f_y,x).doit() #yで偏微分してからxで偏微分

#ヘッセ行列
hesse = Matrix([
[f_xx,f_xy],
[f_yx,f_yy]
])

#ヘッセ行列の行列式
det_hesse = hesse.det()


行列式を出力すると以下のようになりました。

24*x*(4*x + 1) - (8*y - 10)**2

それではこの行列式を用いて極値の判定を行います。


極値の判定

極値の判定は以下の条件に従います。

  • detH(a,b) > 0 かつ f_{xx}(a,b) > 0 ならば ff(a,b) で極小
  • detH(a,b) > 0 かつ f_{xx}(a,b) < 0 ならば ff(a,b) で極大
  • detH(a,b) < 0 ならば ff(a,b)極値を取らない((a,b)は鞍点)
  • detH(a,b) = 0 のときは、個別に調べる必要がある



判定をするコードは以下のようになります。

#判定
for i in range(len(stationary_points)):
    sign_det_hessian = det_hessian.subs([(x,stationary_points[i][0]),(y,stationary_points[i][1])]).evalf()
    if sign_det_hessian > 0:
        sign_f_xx = f_xx.subs([(x,stationary_points[i][0]),(y,stationary_points[i][1])]).evalf()
        if sign_f_xx > 0:
            print('停留点:{0}はf(x,y)の極小点'.format([stationary_points[i][0],stationary_points[i][1]]))
        elif sign_f_xx < 0:
            print('停留点:{0}はf(x,y)の極大点'.format([stationary_points[i][0],stationary_points[i][1]]))
        else:
            print('この方法では判定不可能')
    elif sign_det_hessian < 0:
        print('停留点:{0}はf(x,y)の鞍点'.format([stationary_points[i][0],stationary_points[i][1]]))


結果は以下のようになりました。

停留点:[0, 0]はf(x,y)の鞍点
停留点:[1, 1]はf(x,y)の極小点
停留点:[-3/4 - sqrt(6)/12, -sqrt(6)/8 + 2]はf(x,y)の極大点
停留点:[-3/4 + sqrt(6)/12, sqrt(6)/8 + 2]はf(x,y)の鞍点

まとめ

Sympyを使って比較的簡単に2変数関数の停留点の分類を行うことができました。そんなに数学に詳しくないのでもし誤りがあれば指摘お願いします。
一応すべてのソースコードも載せておきます。
gist.github.com