--- title: "Pythonの最適化ライブラリscipy.optimizeの使い方｜minimize・least_squares・milpと外部ソルバーの選び方【SciPy 1.18】" url: "https://www.issoh.co.jp/tech/details/2049/" published: 2024-02-20 updated: 2026-10-12 categories: ["Python"] publisher: "株式会社一創" --- # Pythonの最適化ライブラリscipy.optimizeの使い方｜minimize・least\_squares・milpと外部ソルバーの選び方【SciPy 1.18】 Pythonで最適化問題を解くとき、ライブラリは「何を最小化したいか」と「変数が連続か整数か」でほぼ決まります。連続な変数の関数最小化やデータへの当てはめならSciPyのscipy.optimizeで完結し、整数変数と業務ルールが増えればPuLPやOR-Toolsのようなモデリング層へ移ります。この記事は、その振り分けの基準を示したうえで、scipy.optimizeの`minimize`・`least_squares`・`linprog`・`milp`を問題の型ごとに使い分ける方法を、実行したコードで解説します。掲載コードと出力は macOS / Python 3.13.14 / SciPy 1.18.1 / NumPy 2.5.4 で実行して取得しました。 ## まとめ：問題の型に応じた関数・ライブラリの選択 | 問題の型 | 使う関数・ライブラリ | 要点 | | ------------- | ------------------------ | -------------- | | 連続変数の関数を最小化 | scipy.optimize.minimize | 既定はBFGS | | 変数に上下限がある | minimize | 既定はL-BFGS-B | | 等式・不等式の制約がある | minimize | 既定はSLSQP | | 残差の二乗和を最小化 | least\_squares | lossで外れ値に強くできる | | モデル式をデータに当てはめ | curve\_fit | 共分散pcovを返す | | 谷が複数ある関数 | differential\_evolution | 初期値でなく範囲を渡す | | 目的関数も制約も線形 | linprog | HiGHSで解く | | 線形で変数が整数 | milp | 小規模ならSciPyで完結 | | 整数変数と業務ルールが多い | PuLP・OR-Tools・Python-MIP | 変数名と式でモデルを書く | | 目的が複数ある | pymoo | パレート解を求める | | 機械学習のパラメータ探索 | Optuna | 試行の打ち切りまで管理 | 表の上8行はSciPy 1.18だけで解けます。外部ライブラリを検討する目安は、制約の添字管理が複雑になった場合と、複数の目的のパレート解の集合を求めたい場合です。以降、SciPyの各関数の使い方を先に示し、最後にSciPyの外へ出す判断を扱います。決定変数・目的関数・制約条件という三要素の整理から始めたい場合は[決定変数・目的関数・制約条件とは？違いと定式化の手順を例題で解説](/tech/details/8321/)が前提になります。 ## 問題の型別のPython最適化ライブラリ選定 ### SciPyによる連続最適化とデータの当てはめ scipy.optimizeは汎用の1関数ではなく、問題の型ごとの関数の集合です。非線形の関数最小化、非線形最小二乗法、大域最適化、線形計画、混合整数線形計画まで揃っており、変数が数十個程度の連続最適化なら追加の依存なしで済みます。線形計画と整数計画のバックエンドはオープンソースのHiGHSで、`linprog`と`milp`の戻り値の`message`にも「HiGHS Status」と出ます。SciPyそのものの位置づけやNumPyとの役割分担は[SciPyとは？読み方・使い方・NumPyとの違いをPythonで解説](/tech/details/3137/)で扱っています。 ### 線形計画・整数計画のモデリングはPuLP・Python-MIP・OR-Tools SciPyの`linprog`と`milp`は、制約を係数行列で渡す形式です。「工場Aの稼働時間」のように変数に名前を付けて式で書きたい規模になったら、モデリング層を持つライブラリに移ります。PuLPは[PyPI](https://pypi.org/project/PuLP/)の説明どおりLPとMILPのモデラーで、CBC・HiGHS・SCIP・Gurobiなど複数のソルバーを呼び分けられます。Python-MIP（PyPI上の名前は`mip`、最新は2026年8月18日公開の2.0.0）はCBCとGurobiとの連携を前提に作られています。スケジューリングのように論理条件が込み入る問題ではOR-ToolsのCP-SATが候補です。 PuLPは2026年9月25日公開の4.0.0で内部がRust実装に置き換わり、互換性が崩れました。検証環境で確認した変化は2つです。 ``` $ pip install pulp==4.0.0 >>> import pulp >>> pulp.listSolvers(onlyAvailable=True) [] >>> pulp.LpVariable("y", 0, 4) TypeError: LpVariable.__init__() takes 2 positional arguments but 4 were given $ pip install "pulp[cbc]==4.0.0" >>> import pulp >>> pulp.listSolvers(onlyAvailable=True) ['COIN_CMD'] ``` `pip install pulp`だけでは使えるソルバーが1つも入らず、CBCは`pulp[cbc]`の追加指定で入ります。変数も`LpVariable(...)`の直接生成ではなく`prob.add_variable("x", 0, 4, cat="Integer")`で作る形に変わりました。3系を前提にした解説記事のコードは4.0.0ではそのまま動かないので、既存環境は版を固定してから移行してください。PuLP 3系までの書き方は[PuLPとは？Pythonで最小化問題を解く手順と3.3系で非推奨になったAPI](/tech/details/3036/)、OR-Toolsは[Google OR-Toolsとは？できること・インストール・実装例をわかりやすく解説](/tech/details/3477/)、商用ソルバーを含む費用とライセンスの比較は[ORツールとは？主要6ツールの比較と選び方【ライセンス・費用・得意問題】](/tech/details/3034/)にまとめています。 ### 凸最適化・多目的最適化・ハイパーパラメータ探索の専用ライブラリ 問題が凸だと分かっているならCVXPYが使えます。数式に近い書き方でモデルを組み、ClarabelやOSQP、HiGHSなどのソルバーへ渡す、凸最適化を中心に混合整数凸計画や幾何計画も扱えるモデリング言語です。目的が「コスト最小」と「品質最大」のように複数あって1つの最適解が決まらない場合は、NSGA-IIなどでパレート解の集合を求めるpymooの領域になります。機械学習モデルのハイパーパラメータ探索はOptunaが定番で、見込みの薄い試行を途中で打ち切る枝刈りの仕組みを持ちます。この系統は勾配法とは設計思想が違うので、[ハイパーパラメータ探索とは？手法の選び方と打ち切り設計を実装目線で解説](/tech/details/16351/)を参照してください。 ## minimizeの最小構成と結果オブジェクトの読み方 ### 目的関数と初期値だけで動く最小コード `minimize`の必須引数は目的関数と初期値`x0`の2つです。目的関数は配列を受け取ってスカラーを返す形にします。 ``` from scipy.optimize import minimize def objective(x): return (x[0] - 3) ** 2 + (x[1] + 1) ** 2 result = minimize(objective, [0.0, 0.0]) print("x =", result.x) print("fun =", result.fun) print("success =", result.success) print("nit =", result.nit, "nfev =", result.nfev) ``` ``` x = [ 3.00000004 -1.00000007] fun = 7.036731787282697e-15 success = True nit = 2 nfev = 9 ``` 戻り値の`OptimizeResult`は、`x`が解、`fun`が目的関数値、`nit`が反復回数、`nfev`が目的関数の評価回数です。解は厳密な3と-1ではなく小数第8位に誤差が残るため、業務ロジックで等値比較をすると破綻します。配列操作の基礎は[NumPyとは？読み方・インストール・基本的な使い方をPythonで解説](/tech/details/3124/)にまとめています。 ### methodを省略したときに選ばれるアルゴリズム [SciPyの公式リファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)は、`method`省略時の規則を「If not given, chosen to be one of BFGS, L-BFGS-B, SLSQP, depending on whether or not the problem has constraints or bounds.」と定めています。制約も範囲もなければBFGS、`bounds`を渡せばL-BFGS-B、`constraints`を渡せばSLSQPです。両方を同時に渡すと検証環境ではSLSQPが選ばれました。`bounds`を足すと実行可能領域が変わるうえ、`method`を省略していればアルゴリズムも入れ替わります。アルゴリズムを固定したいなら`method`を明示してください。 ### successを確認しないままresult.xを使う失敗 `minimize`は収束しなくても例外を投げません。反復上限に達しただけでも`result.x`に値が入り、そのまま下流に流れます。本番コードでは`success`を見て、偽なら`message`とともに止めます。 ``` if not result.success: raise RuntimeError(f"最適化に失敗しました: {result.message}") ``` ただし`success`が真でも、意味は「収束条件を満たした」であって「大域最適に到達した」ではありません。この違いは後半の局所解の節で実例を示します。 ## minimizeのmethodを選ぶ基準とRosenbrock関数での実測 ### 5手法の関数評価回数と到達精度 細長い谷を持ち、素朴な降下法が進みにくいRosenbrock関数で、勾配を渡さない場合の5手法を比べます。 ``` import numpy as np from scipy.optimize import minimize def rosen(x): return (1 - x[0]) ** 2 + 100 * (x[1] - x[0] ** 2) ** 2 x0 = [-1.2, 1.0] for m in ["Nelder-Mead", "Powell", "CG", "BFGS", "L-BFGS-B"]: r = minimize(rosen, x0, method=m) print(f"{m:12s} fun={r.fun:.3e} nfev={r.nfev:4d} nit={r.nit:4d}") ``` ``` Nelder-Mead fun=8.178e-10 nfev= 159 nit= 85 Powell fun=1.787e-26 nfev= 607 nit= 23 CG fun=1.038e-11 nfev= 280 nit= 37 BFGS fun=2.049e-11 nfev= 117 nit= 32 L-BFGS-B fun=9.686e-12 nfev= 138 nit= 37 ``` Powellは桁違いに高精度な代わりに目的関数を607回呼び、BFGSは117回で1e-11台に届いています。目的関数1回がシミュレーション1本に相当する重い問題では、到達精度の桁より`nfev`の差がそのまま実行時間の差になります。 ### Rosenbrock関数での解析勾配jacによる評価回数の削減 BFGS・CG・L-BFGS-Bは勾配を使う手法で、`jac`を渡さないと差分で数値近似します。同じBFGSに解析的な勾配を渡して比べます。 ``` import numpy as np from scipy.optimize import minimize def rosen(x): return (1 - x[0]) ** 2 + 100 * (x[1] - x[0] ** 2) ** 2 def rosen_grad(x): return np.array([ -2 * (1 - x[0]) - 400 * x[0] * (x[1] - x[0] ** 2), 200 * (x[1] - x[0] ** 2), ]) x0 = [-1.2, 1.0] r1 = minimize(rosen, x0, method="BFGS") r2 = minimize(rosen, x0, method="BFGS", jac=rosen_grad) print(f"数値微分 fun={r1.fun:.3e} nfev={r1.nfev} njev={r1.njev}") print(f"解析勾配 fun={r2.fun:.3e} nfev={r2.nfev} njev={r2.njev}") ``` ``` 数値微分 fun=2.049e-11 nfev=117 njev=39 解析勾配 fun=2.535e-15 nfev=39 njev=39 ``` 勾配の評価回数はどちらも39回ですが、数値微分では勾配1回ごとに目的関数を追加で呼ぶため`nfev`が117回に膨らみます。解析勾配なら39回で済み、精度も4桁上がりました。2変数でこの差なので、変数がn個なら差分近似の負担は勾配1回につき目的関数n回前後へ増えます。勾配が書けるなら必ず`jac`に渡してください。 ### 微分できない関数と制約付き問題で選べる手法 目的関数が微分不可能だったり外部プロセスの戻り値でノイズを含んだりするなら、勾配を使わないNelder-MeadかPowellを選びます。等式・不等式の制約を受け付けるのはSLSQP・trust-constr・COBYLA・COBYQAの4手法で、勾配を渡せるならSLSQPかtrust-constr、微分不可能ならCOBYLAかCOBYQAです。dogleg・trust-ncg・trust-exact・trust-krylovは二階微分の`hess`（または`hessp`）を指定しないと`ValueError`で止まります。このうちtrust-ncgとtrust-krylovは`hess='2-point'`の数値近似で動きましたが、doglegとtrust-exactは数値近似を指定してもヘシアン行列を要求して止まりました。同じニュートン法系でもNewton-CGは`jac`だけで動きます。 ## least\_squaresによる非線形最小二乗問題の解き方 ### 残差関数の書き方と戻り値のcost・fun・jac 観測データにモデル式を当てはめる問題は、残差の二乗和を自分で組んで`minimize`に渡すこともできます。しかし`minimize`には、残差のヤコビ行列の構造や外れ値対策の`loss`を直接指定する仕組みがありません。目的関数が残差の二乗和の形をしているなら`least_squares`を使います。渡す関数はスカラーではなく**残差の配列**を返します。指数減衰モデルに、ノイズを加えた40点を当てはめます。 ``` import numpy as np from scipy.optimize import least_squares rng = np.random.default_rng(0) t = np.linspace(0, 10, 40) y = 2.5 * np.exp(-0.4 * t) + 0.3 + rng.normal(0, 0.02, t.size) def residual(p, t, y): a, b, c = p return a * np.exp(-b * t) + c - y r = least_squares(residual, [1.0, 1.0, 0.0], args=(t, y)) print("x =", np.round(r.x, 4)) print("cost =", round(r.cost, 6), " 0.5*sum(fun**2) =", round(0.5 * np.sum(r.fun ** 2), 6)) print("fun.shape =", r.fun.shape, " jac.shape =", r.jac.shape) print("status =", r.status, " nfev =", r.nfev, " optimality =", f"{r.optimality:.1e}") ``` ``` x = [2.5038 0.4066 0.3069] cost = 0.004599 0.5*sum(fun**2) = 0.004599 fun.shape = (40,) jac.shape = (40, 3) status = 2 nfev = 7 optimality = 1.7e-07 ``` 真値a=2.5、b=0.4、c=0.3に近い値に到達し、`nfev`は7でした。`nfev`は数値微分のための呼び出しを含まない数で、残差関数にカウンターを付けて数えると実際には28回呼ばれていました。[least\_squaresの公式リファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least%5Fsquares.html)どおり、既定の`loss='linear'`では`cost`は残差の二乗和の**半分**で、出力の検算でも一致しています。二乗和そのものと取り違えると、モデル比較や残差分散の計算で値が2倍ずれます。`fun`は解での残差ベクトル（40点分）、`jac`は40行3列のヤコビ行列です。`status`は停止理由で、1が`gtol`、2が`ftol`、3が`xtol`、4が`ftol`と`xtol`の両方による収束、0は関数評価の上限到達を表します。`success`は`status`が正のときだけ真になります。 ### lossとf\_scaleによる外れ値の影響の抑制 同じデータの40点中3点に+1.5の外れ値を混ぜ、既定の`loss='linear'`と頑健損失の`soft_l1`を比べます。 ``` import numpy as np from scipy.optimize import least_squares rng = np.random.default_rng(0) t = np.linspace(0, 10, 40) y = 2.5 * np.exp(-0.4 * t) + 0.3 + rng.normal(0, 0.02, t.size) y[[7, 21, 33]] += 1.5 def residual(p, t, y): a, b, c = p return a * np.exp(-b * t) + c - y for loss in ["linear", "soft_l1"]: r = least_squares(residual, [1.0, 1.0, 0.0], args=(t, y), loss=loss, f_scale=0.1) a, b, c = r.x print(f"{loss:8s} a={a:.4f} b={b:.4f} c={c:.4f} nfev={r.nfev}") ``` ``` linear a=2.4313 b=0.3759 c=0.3944 nfev=10 soft_l1 a=2.4986 b=0.4052 c=0.3130 nfev=13 ``` `linear`ではベースラインcが0.3944と真値0.3から31パーセント上振れし、外れ値3点に引っ張られています。`soft_l1`ではcが0.3130まで戻り、`nfev`は10から13へ、数値微分を含む残差関数の呼び出しは40回から46回へ増えただけでした。これはシード0での1回の実測で、外れ値の位置や大きさが変われば差も変わります。`f_scale`は公式リファレンスで「Value of soft margin between inlier and outlier residuals, default is 1.0.」と定義され、残差をどの尺度から抑えるかを決める値です。soft\_l1は境界で残差をはっきり切り分けるのではなく、残差が大きいほど影響を滑らかに弱めます。残差1.5に対する勾配の重みは、`f_scale=1.0`でも二乗損失の約0.55倍に抑えられ、0.1にすると約0.07倍まで下がります。典型的な残差の大きさ（この例ではノイズの標準偏差0.02）を目安に設定してください。文字列で選べる`loss`はlinear・soft\_l1・huber・cauchy・arctanの5つで、自作の関数も渡せます。 ### boundsによる範囲指定とactive\_maskによる境界の確認 パラメータに物理的な範囲があるなら`bounds`で指定します。ここではcの上限を真値より低い0.25に絞り、範囲に張り付く様子を見ます。 ``` import numpy as np from scipy.optimize import least_squares rng = np.random.default_rng(0) t = np.linspace(0, 10, 40) y = 2.5 * np.exp(-0.4 * t) + 0.3 + rng.normal(0, 0.02, t.size) def residual(p, t, y): a, b, c = p return a * np.exp(-b * t) + c - y r = least_squares(residual, [1.0, 1.0, 0.0], args=(t, y), bounds=([0, 0, 0], [10, 2, 0.25])) print("x =", np.round(r.x, 4)) print("active_mask =", r.active_mask) ``` ``` x = [2.5248 0.377 0.25 ] active_mask = [0 0 1] ``` `active_mask`は変数ごとに0なら範囲の内側、-1なら下限、1なら上限に張り付いていることを示します。cが上限0.25に張り付き、その分をaとbが補って真値からずれました。範囲に張り付くこと自体は正しい制約付きの解でも起こりますが、この例のように範囲の設定が真値を外している場合やモデルが合っていない場合にも起こるので、当てはめ後に`active_mask`を確認して見直しの手掛かりにします。 `method`は既定の`trf`のほかに`dogbox`と`lm`があります。`lm`はMINPACKのLevenberg-Marquardt法で、公式リファレンスの説明は「`Doesn't handle bounds and sparse Jacobians.`」です。`bounds`と同時に渡すと`ValueError: Method 'lm' doesn't support bounds.`で止まり、残差の本数が変数の個数以上であることも必要です。範囲付きや疎なヤコビ行列を扱うなら`trf`のままにします。 ### ヤコビ行列によるパラメータ標準誤差の近似 `least_squares`の戻り値には共分散行列が含まれません。推定値の誤差が欲しいときは、解でのヤコビ行列Jから共分散を「Jの転置とJの積の逆行列」に残差分散を掛けて近似します。`curve_fit`の`pcov`と突き合わせます。 ``` import numpy as np from scipy.optimize import least_squares, curve_fit rng = np.random.default_rng(0) t = np.linspace(0, 10, 40) y = 2.5 * np.exp(-0.4 * t) + 0.3 + rng.normal(0, 0.02, t.size) def model(t, a, b, c): return a * np.exp(-b * t) + c r = least_squares(lambda p: model(t, *p) - y, [1.0, 1.0, 0.0]) J = r.jac dof = t.size - r.x.size s2 = 2 * r.cost / dof cov = np.linalg.inv(J.T @ J) * s2 print("least_squares 標準誤差:", np.round(np.sqrt(np.diag(cov)), 5)) popt, pcov = curve_fit(model, t, y, p0=[1.0, 1.0, 0.0]) print("curve_fit 標準誤差:", np.round(np.sqrt(np.diag(pcov)), 5)) ``` ``` least_squares 標準誤差: [0.00951 0.00392 0.00595] curve_fit 標準誤差: [0.00951 0.00392 0.00595] ``` 残差分散は二乗和を自由度（40点−3パラメータ）で割った値で、`cost`が二乗和の半分なので2倍して使っています。5桁まで`curve_fit`と一致しました。この近似は、解の近くでモデルが線形とみなせ、誤差が無相関で等分散、自由度が正、ヤコビ行列の列がフルランクであることが前提で、`loss`に頑健損失を指定したときや解が範囲に張り付いたときは成り立ちません。 ### curve\_fitとleast\_squaresの使い分け [curve\_fitの公式リファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve%5Ffit.html)によれば、既定のmethodは「`'lm' for unconstrained problems and 'trf' if bounds are provided`」です。`curve_fit`は、`lm`では`leastsq`、`trf`と`dogbox`では`least_squares`を内部で呼ぶラッパーで、モデル関数を直接渡せて`pcov`も返します。標準的な当てはめなら`curve_fit`から始めます。`sigma`による重み付けも、`method='trf'`と`loss`の併用も`curve_fit`のままで指定でき、本文の外れ値データに`loss='soft_l1', f_scale=0.1`を渡すと`least_squares`と同じa=2.4986、b=0.4052、c=0.3130が返りました。残差の定義そのものを自分で組みたい場合や、疎なヤコビ行列を指定したい大規模問題で`least_squares`へ降りるのが手戻りの少ない順序です。線形モデルの制約なし最小二乗問題なら、反復型の最適化ソルバーを使わずに`numpy.linalg.lstsq`で解が出ます。 ## 制約付き最適化と、初期値で答えが変わる局所解の対策 ### 等式・不等式制約をSLSQPで解く輸送コスト配分 3拠点から3市場へ配送し、各市場の需要をちょうど満たし、各拠点の供給上限を超えない配分で総コストを最小化します。 ``` import numpy as np from scipy.optimize import minimize cost = np.array([[0.8, 1.2, 1.5], [1.1, 0.7, 1.3], [1.4, 1.0, 0.6]]) demand = [120.0, 80.0, 100.0] supply = [150.0, 100.0, 80.0] def total_cost(x): return float(cost.ravel() @ x) cons = [{"type": "eq", "fun": lambda x, j=j: x.reshape(3, 3)[:, j].sum() - demand[j]} for j in range(3)] cons += [{"type": "ineq", "fun": lambda x, i=i: supply[i] - x.reshape(3, 3)[i].sum()} for i in range(3)] r = minimize(total_cost, np.full(9, 30.0), method="SLSQP", bounds=[(0, None)] * 9, constraints=cons) if not r.success: raise RuntimeError(f"最適化に失敗しました: {r.message}") print("total cost =", round(r.fun, 4)) print(np.round(r.x.reshape(3, 3), 3)) ``` ``` total cost = 226.0 [[120. 0. 0.] [ 0. 80. 20.] [ 0. 0. 80.]] ``` 制約の書式で事故が起きやすいのは2点です。`eq`は「戻り値が0になる」、`ineq`は「戻り値が0以上になる」という意味なので、上限制約は`supply[i] - 使用量`の向きで書きます。もう1つはループ内のラムダ式で、`j=j`のように既定引数で固定しないと全ての制約が最後の`j`を参照します。 ### 初期値を変えると別の解に落ちる問題への対処 勾配法は出発点から近い谷に落ちます。谷が複数ある関数で初期値を3通りに変えます。 ``` import numpy as np from scipy.optimize import minimize, differential_evolution def f(x): return np.sin(3 * x[0]) + 0.1 * (x[0] - 2.0) ** 2 for x0 in [-4.0, 0.0, 4.0]: r = minimize(f, [x0], method="BFGS") print(f"BFGS x0={x0:5.1f} -> x={r.x[0]:7.4f} fun={r.fun:7.4f} success={r.success}") g = differential_evolution(f, bounds=[(-5.0, 5.0)], rng=0) print(f"differential_evolution -> x={g.x[0]:7.4f} fun={g.fun:7.4f} nfev={g.nfev}") ``` ``` BFGS x0= -4.0 -> x=-4.5615 fun= 3.4060 success=True BFGS x0= 0.0 -> x=-0.4685 fun=-0.3770 success=True BFGS x0= 4.0 -> x= 3.6289 fun=-0.7287 success=True differential_evolution -> x= 1.5801 fun=-0.9820 nfev=156 ``` 3つとも`success=True`なのに、目的関数値は3.4060・-0.3770・-0.7287とばらばらで、最良の-0.9820にはどれも届いていません。`differential_evolution`は初期値ではなく探索範囲`bounds`を受け取り、この実行では156回の評価で最小値に到達しました。目的関数の形が読めない問題で、単一の初期値から`minimize`を1回だけ呼ぶ設計は避けてください。各変数の範囲を決められないなら、複数の初期値から回して最良を採る多点スタートにします。範囲を渡す大域最適化には`shgo`と`dual_annealing`もあり、変数が1つで局所的な最小値を探すだけなら`minimize_scalar`も使えますが、多峰性の関数で大域最小値は保証しません。`rng`引数は1.15.0で導入され、それ以前の記事のコードが`seed=`で書かれているのはこのためです。進化計算そのものを使う選択肢は[遺伝的アルゴリズムとは？仕組みとDEAPでのPython実装をわかりやすく解説](/tech/details/3556/)で扱っています。 ## linprogによる線形計画とmilpによる混合整数線形計画 ### 輸送コスト問題でのSLSQPとlinprogの実行時間比較 先ほどの輸送コスト問題は目的関数も制約も線形なので、本来は`linprog`の担当です。同じ問題を初期値なしで解きます。 ``` import numpy as np from scipy.optimize import linprog cost = np.array([[0.8, 1.2, 1.5], [1.1, 0.7, 1.3], [1.4, 1.0, 0.6]]) demand = [120.0, 80.0, 100.0] supply = [150.0, 100.0, 80.0] A_eq = np.zeros((3, 9)) A_ub = np.zeros((3, 9)) for k in range(3): A_eq[k, k::3] = 1.0 A_ub[k, k * 3:(k + 1) * 3] = 1.0 r = linprog(cost.ravel(), A_ub=A_ub, b_ub=supply, A_eq=A_eq, b_eq=demand, bounds=(0, None)) if not r.success: raise RuntimeError(f"最適化に失敗しました: {r.message}") print(r.message) print("total cost =", round(r.fun, 4)) print(np.round(r.x.reshape(3, 3), 3)) ``` ``` Optimization terminated successfully. (HiGHS Status 7: Optimal) total cost = 226.0 [[120. 0. 0.] [ 0. 80. 20.] [ 0. 0. 80.]] ``` 総コストも配分もSLSQPと一致しました。違いは速度と保証です。`timeit`で20回平均を測るコードがこちらです。 ``` import timeit setup = """ import numpy as np from scipy.optimize import minimize, linprog cost = np.array([[0.8, 1.2, 1.5], [1.1, 0.7, 1.3], [1.4, 1.0, 0.6]]) demand = [120.0, 80.0, 100.0] supply = [150.0, 100.0, 80.0] def total_cost(x): return float(cost.ravel() @ x) cons = [{"type": "eq", "fun": lambda x, j=j: x.reshape(3, 3)[:, j].sum() - demand[j]} for j in range(3)] cons += [{"type": "ineq", "fun": lambda x, i=i: supply[i] - x.reshape(3, 3)[i].sum()} for i in range(3)] A_eq = np.zeros((3, 9)); A_ub = np.zeros((3, 9)) for k in range(3): A_eq[k, k::3] = 1.0 A_ub[k, k * 3:(k + 1) * 3] = 1.0 """ t_slsqp = timeit.timeit( 'minimize(total_cost, np.full(9, 30.0), method="SLSQP", bounds=[(0, None)] * 9, constraints=cons)', setup=setup, number=20) / 20 t_lp = timeit.timeit( 'linprog(cost.ravel(), A_ub=A_ub, b_ub=supply, A_eq=A_eq, b_eq=demand, bounds=(0, None))', setup=setup, number=20) / 20 print(f"SLSQP {t_slsqp * 1000:6.1f} ms") print(f"linprog {t_lp * 1000:6.1f} ms") print(f"ratio {t_slsqp / t_lp:.0f}x") ``` 検証環境で18回繰り返すと、SLSQPは82.5〜102.9ミリ秒、`linprog`は4.9〜6.9ミリ秒で、比は13倍から18倍でした。実行時間も速度比も環境と問題の規模で変わるので、この例で`linprog`が1桁速かったという結果として読んでください。加えて`linprog`は線形計画としての最適性が保証され、初期値依存のリスクがありません。 ### milpの整数解と連続解の切り捨ての比較 「トラック2.7台」が許されない問題では、連続解を丸めても最適になる保証はありません。利益5x+4yを最大化する小さな例で、連続解の切り捨てと`milp`を比べます。 ``` import numpy as np from scipy.optimize import linprog, milp, LinearConstraint, Bounds c = np.array([-5.0, -4.0]) # 利益 5x + 4y を最大化したいので符号を反転 A = np.array([[6.0, 4.0], [1.0, 2.0]]) b = np.array([24.0, 6.0]) lp = linprog(c, A_ub=A, b_ub=b, bounds=(0, None)) print("連続解 x =", np.round(lp.x, 4), "利益 =", round(-lp.fun, 4)) x_floor = np.floor(lp.x) print("切り捨て x =", x_floor, "利益 =", float(-c @ x_floor)) ip = milp(c=c, constraints=LinearConstraint(A, -np.inf, b), integrality=np.ones(2), bounds=Bounds(0, np.inf)) print("milp x =", ip.x + 0.0, "利益 =", -ip.fun) ``` ``` 連続解 x = [3. 1.5] 利益 = 21.0 切り捨て x = [3. 1.] 利益 = 19.0 milp x = [4. 0.] 利益 = 20.0 ``` 連続解(3, 1.5)を切り捨てた(3, 1)は利益19ですが、整数の最適解は(4, 0)の利益20で、丸めでは到達できない点にあります。`milp`の[公式リファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html)では`integrality`の値が0なら連続、1なら整数、2なら半連続変数です。`milp`も`linprog`も最小化しか行わないので、最大化は係数の符号を反転し、結果の`fun`も符号を戻して読みます。 ### SciPyから外部ライブラリへ移す判断 最初に見るのは、ソルバーの性能よりモデルの書きやすさです。`linprog`と`milp`は係数行列を自分で組む形式なので、変数間の対応関係や「この設備を動かしたら段取り時間が発生する」といった条件付きの制約の管理が複雑になると、行列の添字管理が先に破綻します。変数名と式でモデルを書けるPuLPやOR-Toolsへ移すのはその時点です。逆に、変数が数十個までの連続最適化やデータの当てはめなら、依存を増やさずSciPyで完結させるほうが運用は軽くなります。求解時間や必要な機能が要件に届かない場合は、接続先のソルバーの性能も比較します。最適化を業務システムへ組み込み、計画の自動作成まで含めて作り込む場合は、要件に合わせて設計する[フルスクラッチ開発](/service/system/fullscratch/)の領域になります。 ## SciPy 1.16〜1.18で変わった点と旧コードの直し方 ### SciPy 1.18系のPython対応バージョン [SciPy 1.18.0のリリースノート](https://docs.scipy.org/doc/scipy/release/1.18.0-notes.html)は「This release requires Python 3.12-3.14 and NumPy 2.0.0 or greater.」としています。2026年8月21日公開の最新版1.18.1は、[GitHubのリリース](https://github.com/scipy/scipy/releases/tag/v1.18.1)によれば新機能のないバグ修正版で、Python 3.15向けのバイナリがPyPIに加わりました。Python 3.11の環境で入る最新は1.17系なので、`pip install -U scipy`で1.18にならないときはインタプリタの版を先に確認してください。 ### 1.18.0で削除された引数 1.18.0では非推奨期間の満了により、`fmin_l_bfgs_b`の`iprint`と`disp`、`nnls`の`atol`が削除されました。関数は残っているので引数を渡さない呼び出しは動きますが、明示している箇所は更新した時点で止まります。 ``` TypeError: fmin_l_bfgs_b() got an unexpected keyword argument 'disp' ``` 移行先は`minimize(..., method="L-BFGS-B")`で、L-BFGS-Bの`options`に`disp`を渡しても「Unknown solver options: disp」の警告が出るだけなので、進捗表示が要るなら`callback`に表示処理を渡します。 ### 1.16.0で変わったleast\_squaresの既定値 [1.16.0のリリースノート](https://docs.scipy.org/doc/scipy/release/1.16.0-notes.html)で、`least_squares`に各反復で呼ばれる`callback`（trfとdogboxのみ）と、数値微分を並列化する`workers`が加わりました。同じ版で`lm`の`max_nfev`の既定値が`100 * n`に変わり、`x_scale`の既定値も1から`None`へ変わっています（`None`のときtrfとdogboxは1、lmは`'jac'`）。1.15以前と同じ`lm`の呼び出しで評価回数や停止位置が変わった場合は、これらの既定値に加え、`lm`の`nfev`から数値微分の呼び出しが除かれた集計方法の変更も確認してください。 ## よくある質問 ### Pythonで最適化問題を解くライブラリはどれを選べばよいですか？ 変数が連続で、関数の最小化かデータへの当てはめならSciPyのscipy.optimizeです。目的関数と制約が線形で整数変数を含む小規模な問題も`milp`で解けます。変数に名前を付けて式でモデルを書きたい規模になったらPuLPやOR-Tools、凸最適化ならCVXPY、目的が複数あるならpymoo、機械学習のハイパーパラメータ探索ならOptunaです。 ### scipy.optimizeで最大化するにはどうすればよいですか？ SciPyに最大化専用の関数はなく、`minimize`も`linprog`も`milp`も最小化だけを行います。目的関数の符号を反転して最小化し、得られた`fun`の符号を戻して読みます。本文の`milp`の例で係数を`[-5.0, -4.0]`と書き、利益を`-ip.fun`で表示しているのがこの操作です。 ### least\_squaresのcostと残差平方和はどう違いますか？ 既定の`loss='linear'`では、`cost`は残差の二乗和の半分です。本文の例では`cost`が0.004599で、`0.5 * np.sum(r.fun ** 2)`と一致しました。このとき残差平方和は`2 * r.cost`、残差分散はそれを自由度（データ点数−パラメータ数）で割った値です。soft\_l1などの頑健損失を指定した場合の`cost`は損失関数を通した値なので、残差平方和は`np.sum(r.fun ** 2)`で計算します。 ### NumPyだけで最適化はできますか？ NumPyに最適化ソルバーは含まれず、`np.minimize`という関数もありません。NumPyで解けるのは`numpy.linalg.lstsq`による線形最小二乗までで、非線形の当てはめや制約付き最適化にはSciPyが必要です。 ### minimizeが収束しないときは何を確認すればよいですか？ まず`result.message`で、反復上限に達したのか精度が出せずに止まったのかを切り分けます。上限到達なら`options={"maxiter": ...}`を増やし、それでも改善しないなら原因は反復回数ではありません。変数間で桁が数桁違うと勾配の推定が荒くなるため、変数を正規化します。初期値を数点変えて結果が割れるなら局所解の問題なので`differential_evolution`か多点スタートへ、問題が線形なら`linprog`へ切り替えます。`linprog`も反復上限や数値的な問題で止まることはあるので、戻り値の`success`と`message`は同じように確認します。 ## 関連記事 - [SciPyとは？読み方・使い方・NumPyとの違いをPythonで解説](/tech/details/3137/) - [数理最適化とは？3要素・アルゴリズム・機械学習との違いをPythonで解説](/tech/details/2976/) - [PuLPとは？Pythonで最小化問題を解く手順と3.3系で非推奨になったAPI](/tech/details/3036/) - [ORツールとは？主要6ツールの比較と選び方【ライセンス・費用・得意問題】](/tech/details/3034/) - [ハイパーパラメータ探索とは？手法の選び方と打ち切り設計を実装目線で解説](/tech/details/16351/) --- 出典: [Pythonの最適化ライブラリscipy.optimizeの使い方｜minimize・least\_squares・milpと外部ソルバーの選び方【SciPy 1.18】]()（株式会社一創）