ほぼPython

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

Flaskで立てたローカルサーバに他のPCやスマホからアクセスする方法

意外としっかり書いてある記事がなかったので書いておきます。

 

Flaskでローカルサーバを立て、そのサーバへ、同じWifiに接続しているPCやスマホからアクセスする方法です。

 

1.app.run(host='0.0.0.0') とする。

host='0.0.0.0'とすることが必要のようです。

 

2.ローカルIPアドレスを調べる。

コマンドプロンプトを開いて ipconfig/all と入力します。表示された中に  IPv4 アドレス という項目があるはずです。そのアドレスをコピーします。

おそらくこんな形です。→192.XXX.XX.X

 

3.他の端末(サーバを立てたPCと同じWifiに接続している機器に限る)からアクセス

先ほどコピーしたローカルIPアドレスとサーバが使っているポートの番号を打てばアクセスできます。(ポート番号は特に指定していなければ5000となります)

こんな形になります→192.XXX.XX.X:5000

 

これで、自分が立てたローカルサーバに同一Wifiに接続している端末からアクセスすることが出来ます。

 

Pythonで多変数関数に対してニュートン法を適用(テイラー展開からアルゴリズム、コードまで)

f:id:short_4010:20190602215202p:plain


ニュートン法とは、解の探索アルゴリズムのひとつです。

 

この方法では、最適化を行いたい目的の関数に対してテイラー級数展開近似を行い、その近似式の停留点(偏微分がすべて0となる点)を更新後の解とします。

 

できるだけわかりやすく書きたいのでテイラー展開からしっかり書いていきます。(PCで見ることをオススメします)

 

テイラー展開

 

テイラー展開とは基準点 \bar{x} から少し離れた点 x における関数値 f(x) を表現するものです。

 

f(x)微分可能な場合、基準点 \bar{x} とそこから \delta x 離れた点  x = \bar{x} + \delta x において、テイラー展開と呼ばれる以下の関係が成り立ちます。

 

\[f(x) = f(\bar{x} + \delta x) = f(\bar{x}) + \sum_{ i = 1 }^{ n } \frac{ \partial f(\bar{x}) }{ \partial x_i }\delta x_i + \frac{1}{2!}\sum_{ i = 1 }^{ n }\sum_{ j = 1 }^{ n } \frac{ \partial ^2 f(\bar{x}) }{ \partial x_i \partial x_j }\delta x_i \delta x_j  + \cdots \]

 

ここで、xn 個の変数を持つベクトルです。

 

そのため例えば  \frac{ \partial f(\bar{x}) }{ \partial x_i }\delta x_i f x_i偏微分した偏導関数 \bar{x} を代入して求めた偏微分係数と  \delta x_i の積を表しています。

 

つまり、  \frac{ \partial f(\bar{x}) }{ \partial x_i }\delta x_i は何らかの数値であり、変数は  \delta x_i だけです。

 

ただこの表現は数学がそう得意ではない人にとって驚くほどわけがわからないと思います。

 

そこで、2変数の場合について書いてみます。

 

2変数の場合のテイラー展開

 

f(x) x_1 と  x_2 の関数であるとします。この場合、上の式は以下のように表せます。

 

\[f(x) =  f(\bar{x} + \delta x) = f(\bar{x}) + \sum_{ i = 1 }^{ 2 } \frac{ \partial f(\bar{x}) }{ \partial x_i }\delta x_i + \frac{1}{2!}\sum_{ i = 1 }^{ 2 }\sum_{ j = 1 }^{ 2 } \frac{ \partial ^2 f(\bar{x}) }{ \partial x_i \partial x_j }\delta x_i \delta x_j  + \cdots \]

 

つまり、

 

f(x) = f(\bar{x}) + \frac{ \partial f(\bar{x}) }{ \partial x_1 }\delta x_1 +  \frac{ \partial f(\bar{x}) }{ \partial x_2 }\delta x_2 + \frac{1}{2!} \frac{ \partial ^2 f(\bar{x}) }{ \partial x_1 \partial x_1 } \delta x_1 \delta x_1  +  \frac{1}{2!} \frac{ \partial ^2 f(\bar{x}) }{ \partial x_1 \partial x_2 }  \delta x_1 \delta x_2 +  \frac{1}{2!}\frac{ \partial ^2 f(\bar{x}) }{ \partial x_2 \partial x_1 }  \delta x_2 \delta x_1 +  \frac{1}{2!}\frac{ \partial ^2 f(\bar{x}) }{ \partial x_2 \partial x_2 }  \delta x_2 \delta x_2 + \cdots

 

となります。

 

多変数として書いた場合よりは幾分わかりやすいですね。

 

テイラー展開を利用した近似

 

さて、このテイラー展開について・・・の部分を省略した式(つまり、2次項まで展開し、3次項以上は無視)を f_{T2}(\bar{x}, \delta x) とします。

 

これは、基準点 \bar{x} において、そこから \delta x だけ離れた点でのfの近似式となります。

 

 \[f_{T2}(\bar{x}, \delta x) = f(\bar{x}) + \sum_{ i = 1 }^{ n } \frac{ \partial f(\bar{x}) }{ \partial x_i }\delta x_i + \frac{1}{2!}\sum_{ i = 1 }^{ n }\sum_{ j = 1 }^{ n } \frac{ \partial ^2 f(\bar{x}) }{ \partial x_i \partial x_j }\delta x_i \delta x_j \]

 

さらに、この式はヤコビ行列  \nabla f(x_k) とヘッセ行列  \nabla ^2 f(x_k) を用いて以下のように書けるようです。(ヤコビ行列とは一階の偏微分係数を並べたもの、ヘッセ行列とは二階の偏微分係数を並べたたものです)

 

 \[f_{T2}(\bar{x}, \delta x) = f(\bar{x}) +  \nabla f(x)^{T}\delta x + \frac{1}{2} \delta x^{T} \nabla ^2 f(\bar{x}) \delta x \]

 

 

ちなみに上付きのTは転置を表しています。僕は実際に確かめてはいませんが、うまくまとめようとすると転置が必要になるのは想像できますね。

 

 上でも書きましたが、この式において変数は \delta x だけであることには注意が必要です。

 

ニュートン法では、この近似式を用います。

 

いざ、ニュートン法

 

 ニュートン法ではまず、答えに近そうな値を適当に選んで、その値を x^{(k)} とします。

 

すると、f_{T2}(x^{(k)}, \delta x^{(k)}) を定めることが出来ます。

 

この式は基準点 x^{(k)} から \delta x^{(k)} だけ離れた位置での関数の近似式でしたね。今回はここから繰り返し処理をするので x^{(k)}\delta x^{(k)} と書いてますが、 \bar{x}\delta x と同じです。

 

次に、この近似式を微分します。何度も書いている通り、この関数において変数は  \delta x^{(k)} だけなので  \delta x^{(k)}微分を行います。

 

\[ \frac{\partial f_{T2}(x_k, \delta x^{(k)})}{\partial \delta x^{(k)}} = \nabla f(x^{(k)}) + \nabla ^2 f(x^{(k)}) \delta x^{(k)} = 0 \]

 

この方程式を満たす  \delta x^{(k)} を用いて、次の基準点 x^{(k+1)} = x^{(k)} + \delta x^{(k)} とします。

 

そしたら、またその基準点を使って f_{T2}(x^{(k+1)}, \delta x^{(k+1)}) を定め、この近似式を微分し、微分がゼロとなる点を基準として、また近似式定めて微分し・・・を繰り返していきます。

 

これにより元の関数の最小値を求めます。

 

ちなみに、なぜ微分がゼロになる点を次の基準点にするかについての説明も軽くしておきます。微分がゼロになる点は停留点と呼ばれ、極大、極小、あるいは鞍点となります。

 

ヘッセ行列が正定値であるときは停留点では極小値をとります。ニュートン法ではヘッセ行列が正定値であるとして計算を進めていくわけです。目的関数の近似式の極小値をたどっていけば、目的関数の最小値にたどり着けそうなイメージはつきますよね。

 

なお、ヘッセ行列が正定値とならないときは最小値にたどり着くことはできません。その点を改善した方法が準ニュートン法であり、ヘッセ行列の代わりにヘッセ行列の近似行列を用います。(BFGS法などがよく知られている)

 

 さて、微分がゼロとなる点および更新点はヘッセ行列の逆行列および  x^{(k+1)} = x^{(k)}+\delta x^{(k)} を用いて以下のように表すことが出来ます。

 

\[ \delta x^{(k)} =  - [\nabla^2 f(x^{(k)})]^{-1}\nabla f(x^{(k)})\]

 

\[ x^{(k+1)} = x^{(k)} - [\nabla^2 f(x^{(k)})]^{-1}\nabla f(x^{(k)})\]

 

さて、やっと理論的な話が終わりました。ここから Python による実装に入っていきます。
 

Python によるアルゴリズムの実装

今回、最小値を求める関数は以下のように定義します。

  • u = x_1 - 0.8
  • v = x_2 -a_0 - a_1 u^{2}  \sqrt{1-u} + a_2u
  • α = -b_0 + b_1  u^{2}\sqrt{1+u} + b_2 u
  • β = \frac{c_0  v^2  (1-c_1v)} { 1 + c_2 u^{2}}
  • f(x_1,x_2) = α  EXP(-β)
  • ただし a = (0.3,0.6,0.2), b = (5,26,3), c = (40,1,10)


驚くほど複雑な関数ですが、この関数はニュートン法に関する論文で紹介された有名なものだそうです。


まずは必要なモジュールをインポートします。

import numpy as np
import sympy
import math

 

次に、関数を定義します。

#最小値を求める関数
x1 = sympy.Symbol("x1")
x2 = sympy.Symbol("x2")

a = [0.3,0.6,0.2]
b = [5,26,3]
c = [40,1,10]
u = x1 - 0.8
v = x2 -(a[0] + a[1] * (u**2) * ((1-u)**(1/2)) - a[2]*u)
alpha = -b[0] + b[1] * (u**2) * ((1+u)**(1/2)) + b[2]*u
beta = c[0] * (v**2) * (1-c[1]*v) / (1 + c[2] * (u**2))
f = alpha * (math.e)**(-beta)


この関数に対して微分を行います。微分にはsympyのDerivativeを利用しました。

#微分
f_x1 = sympy.Derivative(f,x1,1).doit()
f_x2 = sympy.Derivative(f,x2,1).doit()
f_x1_x1 = sympy.Derivative(f_x1,x1,1).doit()
f_x1_x2 = sympy.Derivative(f_x1,x2,1).doit()
f_x2_x2 = sympy.Derivative(f_x2,x2,1).doit()


微分を利用してヤコビ行列とヘッセ行列を返す関数を定義します。これらの関数はある点  (X1,X2) におけるヤコビ行列、ヘッセ行列を返してくれます。

#ヤコビ行列
def yacobi_func(X1,X2):
    f_x1_val = f_x1.subs([(x1,X1),(x2,X2)])
    f_x2_val = f_x2.subs([(x1,X1),(x2,X2)])
    yacobi = np.array([f_x1_val,f_x2_val],dtype="float")
    return yacobi

#ヘッセ行列
def hesse_func(X1,X2):
    f_x1_x1_val = f_x1_x1.subs([(x1,X1),(x2,X2)])
    f_x2_x2_val = f_x2_x2.subs([(x1,X1),(x2,X2)])
    f_x1_x2_val = f_x1_x2.subs([(x1,X1),(x2,X2)])
    hesse = np.array([[f_x1_x1_val,f_x1_x2_val],[f_x1_x2_val,f_x2_x2_val]],dtype="float")
    return hesse


そうしたら、いよいよニュートン法アルゴリズムです。

#ニュートン法のアルゴリズム
def newton_func(x):
    i=0
    max_loop = 100

    while i<max_loop:

        #更新式
        yacobi = yacobi_func(*x)
        hesse = hesse_func(*x)
        hesse_inv = np.linalg.inv(hesse)
        x_new = x - hesse_inv.dot(yacobi)

        #途中経過を表示
        print("K={}  x1:{:.8g}, x2:{:.8g}, ▽f(x):{}".format(i+1,x_new[0],x_new[1],yacobi))

        #終了判定
        f_subs_x = f.subs([(x1,x[0]),(x2,x[1])])
        f_subs_x_new = f.subs([(x1,x_new[0]),(x2,x_new[1])])
        min_diff = 10**(-7)
        if abs(f_subs_x_new - f_subs_x) < min_diff:
            break

        #更新
        i = i+1
        x = x_new

更新式については上で説明した通りです。終了条件については  f(x^{(k)} の値と  f(x^{(k+1)} の値を比較して、その差が小さければほとんど更新がない(つまり、最小値に到達した)と判断しています。


最後に、関数を実行して終わりです。

init_x_type1 = [0.8,0.3]
init_x_type2 = [1.0,0.5]

print("starting point:{}".format(init_x_type1))
newton_func(init_x_type1)
print("")
print("starting point:{}".format(init_x_type2))
newton_func(init_x_type2)
#出力
starting point:[0.8, 0.3]
K=1  x1:0.74230769, x2:0.31153846, ▽f(x):[3. 0.]
K=2  x1:0.73959064, x2:0.31432012, ▽f(x):[-0.08758677 -0.81159701]
K=3  x1:0.73950548, x2:0.31436009, ▽f(x):[ 0.0022284  -0.00647998]
K=4  x1:0.73950546, x2:0.3143601, ▽f(x):[ 3.49570445e-07 -2.66230200e-06]

starting point:[1.0, 0.5]
K=1  x1:1.0748369, x2:0.61133115, ▽f(x):[1.68668034 9.42337253]
K=2  x1:1.0917284, x2:0.72723767, ▽f(x):[1.75837159 2.86926556]
K=3  x1:1.0950945, x2:0.8398244, ▽f(x):[0.80999346 0.97765037]
K=4  x1:1.0879722, x2:0.93425209, ▽f(x):[0.36459031 0.32634061]
K=5  x1:1.0781845, x2:0.95158077, ▽f(x):[0.14481656 0.0477075 ]
K=6  x1:1.0769134, x2:0.95040386, ▽f(x):[ 0.0155495  -0.00333478]
K=7  x1:1.0768959, x2:0.95040621, ▽f(x):[2.06073682e-04 1.00748357e-05]

この結果からわかるように、初期値により最終的な (x_1,x_2) は異なります。

どちらが正しい最小点かというと、初期値が(0.8, 0.3)の時です。

ニュートン法では、初期値は最小点の付近の値を与えないとうまく収束しません。

最小点付近の値は関数を可視化したり、等高線を表示することで検討を付けます。

ちなみに今回の関数のグラフはこのようになりました。

f:id:short_4010:20190602213416p:plain

たしかに、(0.8, 0.3)付近で最小となっていることがわかりますね。

さいごに

長くなりましたが、これでニュートン法の理論の説明から実装までが終わりとなります。

この記事がどこかの誰かの役に立つことを願っています。

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