SciPyによる最適化(scipy.optimize)②|SciPy入門 #6

f:id:lib-arts:20190225182115g:plain

SciPyについて色々と話題になり面白そうだったので公式チュートリアルを元にまとめています。

SciPy Tutorial — SciPy v1.2.1 Reference Guide
#5ではscipy.optimizeから制約条件のない際の最適化に関して取り扱いました。

 #6では制約条件がある場合の最適化や最小二乗法などについて取り扱っていければと思います。

Optimization (scipy.optimize) — SciPy v1.2.1 Reference Guide
以下目次になります。

1. Constrained minimization of multivariate scalar functions (minimize)
2. Least-squares minimization (least_squares)
3. まとめ

 

1. Constrained minimization of multivariate scalar functions (minimize)
まずは制約条件の設定の仕方について色々とまとまっています。
・定義域の制約
 0 \leq x_{0} \leq 1 -0.5 \leq x_{1} \leq 2.0のような最適化対象の変数の定義域の制約についてはBoundsオブジェクトで表します。

from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])

これによって最適化対象の変数の定義域を指定することができます。

 

・線形式の制約
 x_{0} + 2x_{1} \leq 1 2x_{0} + x_{1} = 1のような線形式は下記のように行列表記に直した上で考えます。

f:id:lib-arts:20190302114853p:plain
等号を表すにあたって、同じ値で挟んで表現していることに着目です。この変換に基づいて、LinearConstraintオブジェクトに値を設定します。

from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

 

非線形の制約

f:id:lib-arts:20190302114747p:plain
上のような制約についてはNonlinearConstraintオブジェクトを用いて下記のように表します。

def cons_f(x):
    return [x[0]**2 + x[1], x[0]**2 - x[1]]
def cons_J(x):
    return [[2*x[0], 1], [2*x[0], -1]]
def cons_H(x, v):
    return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)

 

最適化問題の計算
#5でまとめたminimizeにおいて上記でまとめた制約条件をconstraintsやboundsで与えてやることで、制約条件ありの最適化を行うことができます。

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

 


2. Least-squares minimization (least_squares)
SciPyは制約つき非線形最小二乗問題(constrained nonlinear least-squares problems)を解くこともできると冒頭で言及した上で下記で問題定義を行っています。

f:id:lib-arts:20190302120330p:plain

上記はあくまで一般的な表記なので、実際に具体的な問題で解説されています。

f:id:lib-arts:20190302120901p:plain
具体的な問題の表記は上記なのですが、これは酵素反応(Enzyme Reaction)をモデル化した数式のようです。ここで解きたいフィッティングの問題は、独立変数(independent variable)のuに基づいて測定値(measurement values)を予測するフィッティングの問題です。通常のフィッティングでは独立変数として使われることの多いxがここでは推定するモデルのパラメータになっていることに注意です。この問題の解法はleast_squaresを用いることで下記で解くことができます。

from scipy.optimize import least_squares
def model(x, u):
    return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])
def fun(x, u, y):
    return model(x, u) - y
def jac(x, u, y):
    J = np.empty*1
    den = u ** 2 + x[2] * u + x[3]
    num = u ** 2 + x[1] * u
    J[:, 0] = num / den
    J[:, 1] = x[0] * u / den
    J[:, 2] = -x[0] * num * u / den ** 2
    J[:, 3] = -x[0] * num / den ** 2
    return J
u = np.array([4.0, 2.0, 1.0, 5.0e-1, 2.5e-1, 1.67e-1, 1.25e-1, 1.0e-1,
            8.33e-2, 7.14e-2, 6.25e-2])
y = np.array([1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2,
            4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2])
x0 = np.array([2.5, 3.9, 4.15, 3.9])
res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=1)
print(res.x)

また、チュートリアルのその下のコードを実行することで下記のようにデータとフィッティングされたモデルを図示することが可能です。

f:id:lib-arts:20190302121853p:plain

元々のフィッティングする数式が難しそうに見えますが、uに0を入れた際にyが0なのと、uを限りなく大きくするとx0の値に近づくことを意識するとなんとなくグラフの形にを把握できるかと思います。



3. まとめ
#6では制約ありの最適化問題と最小二乗問題の最適化について取り扱いました。
#7では残りのトピックである1変数の関数の最適化や求根などを取り扱い、最適化(scipy.optimize)のチュートリアルについては一旦一区切りとできればと思います。

 

*1:u.size, x.size