SciPy (サイパイと発音) は、C および Fortran ライブラリも含まれる numpy ベースの数学パッケージです。 SciPy は、対話型の Python セッションを MATLAB、IDL、Octave、R、SciLab などの完全なデータ サイエンス環境に変えます。
この記事では、数理計画法の基本テクニック、つまり scipy.optimize パッケージを使用して複数の変数のスカラー関数の条件付き最適化問題を解く方法を見ていきます。 制約のない最適化アルゴリズムについては、すでに説明されています。 。 scipy 関数に関するより詳細な最新のヘルプは、help() コマンド、Shift+Tab、または .
導入
scipy.optimize パッケージの条件付きおよび制約なしの両方の最適化問題を解決するための共通インターフェイスは、関数によって提供されます。 minimize()。 ただし、すべての問題を解決するための普遍的な方法はないことが知られているため、適切な方法の選択は常に研究者の肩にかかっています。
適切な最適化アルゴリズムは、関数の引数を使用して指定されます。 minimize(..., method="").
複数の変数の関数の条件付き最適化には、次のメソッドの実装が利用できます。
trust-constr— 信頼領域内の極小値を検索します。 , ;SLSQP— 制約付き逐次二次計画法、ラグランジュ系を解くためのニュートン法。 .TNC- 切り捨てられたニュートン 制約があり、反復回数が制限されており、多数の独立変数を持つ非線形関数に適しています。 .L-BFGS-B— Broyden-Fletcher-Goldfarb-Shanno チームによるメソッド。ヘッセ行列からのベクトルの部分的な読み込みによりメモリ消費量が削減されて実装されています。 , .COBYLA— MARE 線形近似による制約付き最適化、線形近似による制約付き最適化 (勾配計算なし)。 .
選択した方法に応じて、問題を解決するための条件と制限が異なります。
- クラス オブジェクト
Boundsメソッド L-BFGS-B、TNC、SLSQP、trust-constr の場合。 - リスト
(min, max)同じメソッドの場合、L-BFGS-B、TNC、SLSQP、trust-constr。 - オブジェクトまたはオブジェクトのリスト
LinearConstraint,NonlinearConstraintCOBYLA、SLSQP、trust-constr メソッドの場合。 - 辞書または辞書のリスト
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}COBYLA、SLSQP メソッドの場合。
記事の概要:
1) 制約をオブジェクトとして指定して、信頼領域 (method=”trust-constr”) で条件付き最適化アルゴリズムの使用を検討します。 Bounds, LinearConstraint, NonlinearConstraint ;
2) 辞書形式で制限を指定した最小二乗法 (method = "SLSQP") を使用した逐次プログラミングを検討します。 {'type', 'fun', 'jac', 'args'};
3) Web スタジオの例を使用して、製品の最適化の例を分析します。
条件付き最適化メソッド = "trust-constr"
メソッドの実装 trust-constr に基づく 等式の制約を伴う問題および 不等式の形の制約を伴う問題の場合。 どちらの方法も、信頼領域内の極小値を見つけるアルゴリズムによって実装されており、大規模な問題に適しています。
一般形式での最小値を求める問題の数学的定式化:
厳密な等価制約の場合、下限は上限と等しく設定されます。 .
一方向制約の場合、上限または下限が設定されます np.inf 対応する記号が付いています。
XNUMX つの変数の既知のローゼンブロック関数の最小値を見つける必要があるとします。
この場合、その定義ドメインには次の制限が設定されます。
私たちの場合、この時点で独自の解決策があります。 、最初と XNUMX 番目の制限のみが有効です。
制限を下から上に見て、scipy で制限を記述する方法を見てみましょう。
制限 и Bounds オブジェクトを使用して定義しましょう。
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])制限 и 線形形式で書いてみましょう。
これらの制約を LinearConstraint オブジェクトとして定義しましょう。
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])最後に行列形式の非線形制約:
この制約のヤコビアン行列と、ヘッセ行列と任意のベクトルの線形結合を定義します。 :
これで、非線形制約をオブジェクトとして定義できるようになりました。 NonlinearConstraint:
from scipy.optimize import 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]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)サイズが大きい場合は、行列をスパース形式で指定することもできます。
from scipy.sparse import csc_matrix
def cons_H_sparse(x, v):
return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_sparse)またはオブジェクトとして LinearOperator:
from scipy.sparse.linalg import LinearOperator
def cons_H_linear_operator(x, v):
def matvec(p):
return np.array([p[0]*2*(v[0]+v[1]), 0])
return LinearOperator((2, 2), matvec=matvec)
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
jac=cons_J, hess=cons_H_linear_operator)ヘシアン行列を計算する場合 多くの労力が必要なので、クラスを使用できます 。 次の戦略が利用可能です。 BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())ヘッセ行列は、有限差分を使用して計算することもできます。
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')制約のヤコビアン行列は、有限差分を使用して計算することもできます。 ただし、この場合、有限差分を使用してヘッセ行列を計算することはできません。 ヘッセ行列は関数として定義するか、HessianUpdateStrategy クラスを使用して定義する必要があります。
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())最適化問題の解決策は次のようになります。
from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod
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)`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]必要に応じて、ヘッセ行列を計算する関数は LinearOperator クラスを使用して定義できます。
def rosen_hess_linop(x):
def matvec(p):
return rosen_hess_prod(x, p)
return LinearOperator((2, 2), matvec=matvec)
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)または、ヘッセ行列とパラメータを介した任意のベクトルの積 hessp:
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)あるいは、最適化されている関数の一次導関数と二次導関数を近似的に計算することもできます。 たとえば、ヘッセ行列は次の関数を使用して近似できます。 SR1 (準ニュートン近似)。 勾配は有限差分で近似できます。
from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(),
constraints=[linear_constraint, nonlinear_constraint],
options={'verbose': 1}, bounds=bounds)
print(res.x)条件付き最適化手法="SLSQP"
SLSQP メソッドは、次の形式で関数を最小化する問題を解決するように設計されています。
どこ и — 等式または不等式の形式で制限を記述する式のインデックスのセット。 — 関数の定義領域の下限と上限のセット。
線形制約と非線形制約は、キーを含む辞書の形式で記述されます。 type, fun и jac.
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
1 - x [0] ** 2 - x [1],
1 - x [0] ** 2 + x [1]]),
'jac': lambda x: np.array ([[- 1.0, -2.0],
[-2 * x [0], -1.0],
[-2 * x [0], 1.0]])
}
eq_cons = {'type': 'eq',
'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
'jac': lambda x: np.array ([2.0, 1.0])
}最小値の検索は次のように実行されます。
x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
bounds=bounds)
print(res.x)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]最適化の例
第 XNUMX の技術構造への移行に関連して、小規模ながら安定した収入をもたらしてくれる Web スタジオを例に、生産の最適化を見てみましょう。 自分自身が XNUMX 種類の製品を生産する調理室のディレクターであると想像してみましょう。
- x0 - ランディング ページの販売、10 tr から。
- x1 - 企業ウェブサイト、20 tr から。
- x2 - オンラインストア、30 tr から。
私たちのフレンドリーなワーキングチームには、ジュニア XNUMX 名、ミドル XNUMX 名、シニア XNUMX 名が含まれています。 彼らの毎月の労働時間資金:
- ジュネス:
4 * 150 = 600 чел * час, - ミドル:
2 * 150 = 300 чел * час, - 上院議員:
150 чел * час.
最初に参加可能なジュニアは、タイプ (x0、x1、x2)、中級者 - (10、20、30)、上級者 - (7、15、20) の 5 つのサイトの開発と展開に (10、15、XNUMX) 時間を費やします。 ) 人生で最高の時間を過ごす時間。
通常の取締役と同様に、私たちも毎月の利益を最大化したいと考えています。 成功への最初のステップは、目的関数を書き留めることです value XNUMXか月あたりに生産される製品からの収入としては、次のようになります。
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]これはエラーではありません。最大値を検索すると、目的関数は反対の符号で最小化されます。
次のステップは、従業員の過重労働を禁止し、労働時間の制限を導入することです。
同等のもの:
ineq_cons = {'type': 'ineq',
'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
300 - 7 * x [0] - 15 * x [1] - 20 * x[2],
150 - 5 * x [0] - 10 * x [1] - 15 * x[2]])
}正式な制限として、積の出力は正のみでなければならないということです。
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])そして最後に、最も楽観的な仮定は、低価格と高品質のおかげで、満足した顧客が常に列を作っているということです。 制約付き最適化問題を解くことに基づいて、毎月の生産量を自分たちで選択できます。 scipy.optimize:
x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)[7.85714286 5.71428571 3.57142857]大まかに整数に四捨五入して、製品の最適な配分による漕ぎ手の月間負荷を計算してみましょう。 x = (8, 6, 3) :
- ジュネス:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - ミドル:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - 上院議員:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
結論: ディレクターが当然の最大限の成果を得るには、月に 8 つのランディング ページ、6 つの中規模サイト、3 つのストアを作成するのが最適です。 この場合、上級者は機械から顔を上げずに耕す必要があり、中間者の負荷は約2/3、下級者は半分以下になります。
まとめ
この記事では、パッケージを操作するための基本的なテクニックについて概説します。 scipy.optimize、条件付き最小化問題を解くために使用されます。 個人的に私が使っているのは scipy 純粋に学術目的のため、ここで挙げた例は非常に滑稽な性質を持っています。
たとえば、I.L. Akulich の著書「例と問題における数学的プログラミング」には、多くの理論と仮想例が記載されています。 より本格的なアプリケーション scipy.optimize 一連の画像から 3D 構造を構築するには (）でご覧いただけます .
主な情報源は、 このセクションおよび他のセクションの翻訳に貢献したい人 scipy へようこそ .
感謝 出版物の準備に参加するため。
出所： habr.com