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 で実行して取得しました。

まとめ：問題の型に応じた関数・ライブラリの選択

表の上8行はSciPy 1.18だけで解けます。外部ライブラリを検討する目安は、制約の添字管理が複雑になった場合と、複数の目的のパレート解の集合を求めたい場合です。以降、SciPyの各関数の使い方を先に示し、最後にSciPyの外へ出す判断を扱います。決定変数・目的関数・制約条件という三要素の整理から始めたい場合は決定変数・目的関数・制約条件とは？違いと定式化の手順を例題で解説が前提になります。

問題の型別のPython最適化ライブラリ選定

SciPyによる連続最適化とデータの当てはめ

scipy.optimizeは汎用の1関数ではなく、問題の型ごとの関数の集合です。非線形の関数最小化、非線形最小二乗法、大域最適化、線形計画、混合整数線形計画まで揃っており、変数が数十個程度の連続最適化なら追加の依存なしで済みます。線形計画と整数計画のバックエンドはオープンソースのHiGHSで、linprogとmilpの戻り値のmessageにも「HiGHS Status」と出ます。SciPyそのものの位置づけやNumPyとの役割分担はSciPyとは？読み方・使い方・NumPyとの違いをPythonで解説で扱っています。

線形計画・整数計画のモデリングはPuLP・Python-MIP・OR-Tools

SciPyのlinprogとmilpは、制約を係数行列で渡す形式です。「工場Aの稼働時間」のように変数に名前を付けて式で書きたい規模になったら、モデリング層を持つライブラリに移ります。PuLPはPyPIの説明どおり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、OR-ToolsはGoogle OR-Toolsとは？できること・インストール・実装例をわかりやすく解説、商用ソルバーを含む費用とライセンスの比較はORツールとは？主要6ツールの比較と選び方【ライセンス・費用・得意問題】にまとめています。

凸最適化・多目的最適化・ハイパーパラメータ探索の専用ライブラリ

問題が凸だと分かっているならCVXPYが使えます。数式に近い書き方でモデルを組み、ClarabelやOSQP、HiGHSなどのソルバーへ渡す、凸最適化を中心に混合整数凸計画や幾何計画も扱えるモデリング言語です。目的が「コスト最小」と「品質最大」のように複数あって1つの最適解が決まらない場合は、NSGA-IIなどでパレート解の集合を求めるpymooの領域になります。機械学習モデルのハイパーパラメータ探索はOptunaが定番で、見込みの薄い試行を途中で打ち切る枝刈りの仕組みを持ちます。この系統は勾配法とは設計思想が違うので、ハイパーパラメータ探索とは？手法の選び方と打ち切り設計を実装目線で解説を参照してください。

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で解説にまとめています。

methodを省略したときに選ばれるアルゴリズム

SciPyの公式リファレンスは、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の公式リファレンスどおり、既定の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の公式リファレンスによれば、既定の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実装をわかりやすく解説で扱っています。

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の公式リファレンスではintegralityの値が0なら連続、1なら整数、2なら半連続変数です。milpもlinprogも最小化しか行わないので、最大化は係数の符号を反転し、結果のfunも符号を戻して読みます。

SciPyから外部ライブラリへ移す判断

最初に見るのは、ソルバーの性能よりモデルの書きやすさです。linprogとmilpは係数行列を自分で組む形式なので、変数間の対応関係や「この設備を動かしたら段取り時間が発生する」といった条件付きの制約の管理が複雑になると、行列の添字管理が先に破綻します。変数名と式でモデルを書けるPuLPやOR-Toolsへ移すのはその時点です。逆に、変数が数十個までの連続最適化やデータの当てはめなら、依存を増やさずSciPyで完結させるほうが運用は軽くなります。求解時間や必要な機能が要件に届かない場合は、接続先のソルバーの性能も比較します。最適化を業務システムへ組み込み、計画の自動作成まで含めて作り込む場合は、要件に合わせて設計するフルスクラッチ開発の領域になります。

SciPy 1.16〜1.18で変わった点と旧コードの直し方

SciPy 1.18系のPython対応バージョン

SciPy 1.18.0のリリースノートは「This release requires Python 3.12-3.14 and NumPy 2.0.0 or greater.」としています。2026年8月21日公開の最新版1.18.1は、GitHubのリリースによれば新機能のないバグ修正版で、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のリリースノートで、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で解説
• 数理最適化とは？3要素・アルゴリズム・機械学習との違いをPythonで解説
• PuLPとは？Pythonで最小化問題を解く手順と3.3系で非推奨になったAPI
• ORツールとは？主要6ツールの比較と選び方【ライセンス・費用・得意問題】
• ハイパーパラメータ探索とは？手法の選び方と打ち切り設計を実装目線で解説