pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

トップ > Python > Python（SciPy）で学ぶ最適化問題（minimizeの使い方）

Python（SciPy）で学ぶ最適化問題（minimizeの使い方）

Python Python-SciPy Python-SciPy-最適化

数学の具体的な問題にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonの科学技術計算ライブラリであるSciPyを使った最適化問題の例として、関数の極値を求める問題を扱います。また補助的に数値計算やグラフの描画を用いるので、NumPy、matplotlibの使い方の復習にもなるかと思います。



関数の極値：１変数の場合

まず手始めに、次の１変数関数

\ \displaystyle{\qquad f(x)=-x^2+x^4}

が極値（ f'(x)=0）をとる点を求めてみましょう。導関数は

\ \displaystyle{\qquad f'(x)=-2x+4x^3=2x(2x^2-1)}

ですから、関数 f(x) は

\ \displaystyle{\qquad x_0=0 \qquad x_\pm =\pm \frac{1}{\sqrt{2}} }

において極値をとります。グラフを描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt 
%matplotlib inline

x_range = [-1.5, 1.5] 
y_range = [-1, 2]

x = np.linspace(-1.5, 1.5, 60)
y = - x**2 + x**4

plt.plot(x,y)
plt.grid(True)
plt.xlim(x_range)
plt.ylim(y_range)
plt.show()

より

と表示されます。特に、極小（かつ最小）となる点が

\ \displaystyle{\qquad x_\pm =\pm \frac{1}{\sqrt{2}} =\pm 0.7071\dots}

であることがわかります。これらの点での関数の値は

\ \displaystyle{\qquad f(x_\pm) -\frac{1}{\sqrt{2}^2} +\frac{1}{\sqrt{2}^4} =-\frac{1}{4}= -0.25}

です。さて、適当に設定した初期位置の直近にある極小点は、SciPyの組み込み関数 minimize を使って求めることができます。いま初期位置として x=1.0 ととると、直近の極小点は x=x_+=0.7071\dots のはずです。実際

from scipy.optimize import minimize

def f(x):
    y = - x**2 + x**4
    return y

x_init = 1.0

minimize( f, x_init )

とすると

      fun: -0.24999999999984324
 hess_inv: array([[0.24908418]])
      jac: array([-1.08964741e-06])
  message: 'Optimization terminated successfully.'
     nfev: 14
      nit: 5
     njev: 7
   status: 0
  success: True
        x: array([0.7071065])

といった出力を得ます。出力の最後（.x）と最初（.fun）がそれぞれ直近の極小点の位置と、その点での関数の値です。

minimize( f, x_init ).x

minimize( f, x_init ).fun

上で求めていた結果と合っていることがわかります。
　





平面上の３点からの距離和の最小点（フェルマー点）：２変数の場合

問題を少し複雑にしてみましょう。最適化問題の初等的な具体例として、フェルマー点を求める問題を考えてみます。フェルマー点とは、 xy-平面上の適当に与えられた３点（例としてここでは

\ \displaystyle{\qquad {\rm O}=(0,0) \qquad\qquad {\rm P}=(1,0) \qquad\qquad {\rm Q}=(0,1) }

とします）を考え、点 {\rm C}=(a,b) をとるとき、各点O、P、Qから点Cまでの距離の和

\ \displaystyle{\qquad D(a,b)=|{\rm OC}|+|{\rm PC}|+|{\rm QC}| }

を最小にするような点です。フェルマー点についての詳しい解説は

ja.wikipedia.org

など、他の方による解説記事にまかせることにします。関数 D(a,b) の具体形はいまの場合

\ \displaystyle{\qquad D(a,b)=\sqrt{a^2+b^2}+\sqrt{(a-1)^2+b^2}+\sqrt{a^2+(b-1)^2} }

で与えられます。この関数を極小にするような点 (a,b) を求めるのですが、このような場合にもSciPyの minimize を使うことができます。いま初期位置として

\ \displaystyle{\qquad (a_{\rm init},b_{\rm init})=(0.3,\ 0.3) }

などと適当にとって

import numpy as np
from scipy.optimize import minimize

def D(c):
    d = np.sqrt(c[0]**2+c[1]**2) + np.sqrt((c[0]-1)**2+c[1]**2) \
        + np.sqrt(c[0]**2+(c[1]-1)**2)
    return d

c_init = [.3, .3]

minimize(D, c_init)

とすると

\ \displaystyle{\qquad (a_{\rm min},b_{\rm min})=(0.2113\dots,\ 0.2113\dots) }

が求める点であることが出力されます。初等幾何によって示すことのできるフェルマー点の性質から

\ \displaystyle{\qquad \angle{\rm OCP}=\angle{\rm PCQ}=\angle{\rm QCO}=120^\circ }

となることが知られていますが、実際、ベクトル

\ \displaystyle{\qquad \begin{aligned}\overrightarrow{{\rm CO}}&=(-a_{\rm min},-b_{\rm min}) \\ \overrightarrow{{\rm CP}}&=(1-a_{\rm min},-b_{\rm min}) \\ \overrightarrow{{\rm CQ}}&=(-a_{\rm min},1-b_{\rm min})\end{aligned} }

を使って（ a_{\rm min}=b_{\rm min}=:c_{\rm min} として）

\ \displaystyle{\qquad \begin{aligned} \cos\angle{\rm OCP} &=\frac{\overrightarrow{{\rm CO}}\cdot\overrightarrow{{\rm CP}}} {|\overrightarrow{{\rm CO}}|\,|\overrightarrow{{\rm CP}}|} =\frac{2c_{\rm min}^2-c_{\rm min}} {\sqrt{c_{\rm min}^2+c_{\rm min}^2}\sqrt{(1-c_{\rm min})^2+c_{\rm min}^2}} \sim -0.5=\cos 120^\circ \\ \cos\angle{\rm PCQ} &=\frac{\overrightarrow{{\rm CP}}\cdot\overrightarrow{{\rm CQ}}} {|\overrightarrow{{\rm CP}}|\,|\overrightarrow{{\rm CQ}}|} =\frac{2c_{\rm min}^2-2c_{\rm min}} {\sqrt{(1-c_{\rm min})^2+c_{\rm min}^2}\sqrt{c_{\rm min}^2+(1-c_{\rm min})^2}} \sim -0.5=\cos 120^\circ \\ \cos\angle{\rm QCO} &=\frac{\overrightarrow{{\rm CQ}}\cdot\overrightarrow{{\rm CO}}} {|\overrightarrow{{\rm CQ}}|\,|\overrightarrow{{\rm CO}}|} =\frac{2c_{\rm min}^2-c_{\rm min}} {\sqrt{c_{\rm min}^2+(1-c_{\rm min})^2}\sqrt{c_{\rm min}^2+c_{\rm min}^2}} \sim -0.5=\cos 120^\circ \end{aligned} }

から確認することができます。

c_min = 0.211325
cos1 = (2*c_min**2-c_min)/(np.sqrt(2*c_min**2)*np.sqrt((1-c_min)**2+c_min**2))
cos2 = (2*c_min**2-2*c_min)/(np.sqrt((1-c_min)**2+c_min**2)*np.sqrt(c_min**2+(1-c_min)**2))
cos3 = (2*c_min**2-c_min)/(np.sqrt(c_min**2+(1-c_min)**2)*np.sqrt(c_min**2+c_min**2))
[cos1, cos2, cos3]

　



拘束条件のある場合

関数

\ \displaystyle{\qquad g(x,y)= y+xy }

を曲線

\ \displaystyle{\qquad c(x,y)=x^2+y-1=0 }

の上の点に制限して、その極値問題を考えましょう。いまの場合、曲線の式を y について解いてしまって

\ \displaystyle{\qquad y=-x^2+1 }

として関数 g に代入すれば

\ \displaystyle{\qquad \begin{aligned} g(x,y=-x^2+1)&= -x^2+1+x(-x^2+1)\\ &=-x^3-x^2+x+1\\ &= -(x-1)(x+1)^2 \end{aligned} }

であり、極値点を調べるために x について微分すると

\ \displaystyle{\qquad \begin{aligned} g'(x,y=-x^2+1)&= -(x+1)^2-2(x-1)(x+1)\\ &=-(x+1)\left\{(x+1)+2(x-1)\right\}\\ &=-(x+1)(3x-1) \end{aligned} }

から

\ \displaystyle{\qquad(x,y)=(-1,0), \quad (\frac13,\ \frac89) }

で極値をとることがわかります。実際に関数を図示してみましょう。

import matplotlib.pyplot as plt 
%matplotlib inline

x_range = [-3, 3] 
y_range = [-6, 6]

x = np.linspace(-3,3, 50)
y = - (x-1)*(x+1)**2

plt.plot(x,y)
plt.grid(True)
plt.xlim(x_range)
plt.ylim(y_range)
plt.show()

より

といったグラフを得ます。特に極小点は

\ \displaystyle{\qquad(x,y)=(-1,0) }

であることがわかります。いまの場合には拘束条件が簡単だったので式変形できてしまいましたが、より一般的にはラグランジュ未定乗数 \lambda を導入して新たに関数 \widetilde{g} を

\ \displaystyle{\qquad \begin{aligned} \widetilde{g}(x,y;\lambda) &=g(x,y)-\lambda\,c(x,y)\\ &=y+xy -\lambda\left(x^2+y-1\right) \end{aligned} }

と定義し、 x,y,\lambda を独立変数とみなして偏導関数

\ \displaystyle{\qquad \begin{aligned} \frac{\partial \widetilde{g}}{\partial x}&=y-2\lambda x \\ \frac{\partial \widetilde{g}}{\partial y}&=1+x-\lambda \\ \frac{\partial \widetilde{g}}{\partial \lambda}&=-x^2-y+1 \end{aligned} }

を考えて、これらについて連立式

\ \displaystyle{\qquad \begin{aligned} \frac{\partial \widetilde{g}}{\partial x}&=0 \\ \frac{\partial \widetilde{g}}{\partial y}&=0 \\ \frac{\partial \widetilde{g}}{\partial \lambda}&=0 \end{aligned} }

を課す、という方法によって極値点を求めることができます（ラグランジュの未定乗数法）。このような拘束条件付きの場合にも minimize 関数は対応していて

import numpy as np
from scipy.optimize import minimize

def g(x):
    z = x[1] + x[0] * x[1]
    return z

def con(x):
    c = x[0]**2 + x[1] -1
    return c

cons = {'type':'eq', 'fun': con}

x_init = [-1., -1.]

minimize(g, x_init, constraints=cons)

によって確かに

\ \displaystyle{\qquad(x,y)=(-1,0) }

が極小であることが出力されます。







いかがだったでしょうか。今回は、最適化問題をSciPyで扱う方法、特に組み込み関数 minimize の使い方の初歩について解説してみました。


キーワード：Python、SciPy、最適化問題

プライバシーポリシー