SciPy-1-12-中文文档-十-
SciPy 1.12 中文文档(十)
scipy.optimize.OptimizeResult
class scipy.optimize.OptimizeResult
表示优化结果。
注意事项
根据具体使用的求解器,OptimizeResult
可能不包含此处列出的所有属性,并且可能有其他未列出的属性。由于该类本质上是 dict 的子类,带有属性访问器,可以使用 OptimizeResult.keys
方法查看可用的属性。
属性:
xndarray
优化的解。
success布尔值
优化器是否成功退出。
status整数
优化器的终止状态。其值取决于底层求解器。详细信息请参考 message。
message字符串
终止原因的描述。
fun, jac, hess: ndarray
目标函数的值,其雅可比矩阵及海森矩阵的值(如果可用)。这些海森矩阵可能是近似值,请参阅相关函数的文档。
hess_inv对象
目标函数海森矩阵的逆;可能是一个近似值。并非所有求解器都支持。此属性的类型可以是 np.ndarray 或 scipy.sparse.linalg.LinearOperator。
nfev, njev, nhev整数
目标函数及其雅可比矩阵和海森矩阵的评估次数。
nit整数
优化器执行的迭代次数。
maxcv浮点数
最大约束违规。
方法
__getitem__ |
x.getitem(y) <==> x[y] |
---|---|
__len__ (/) |
返回 len(self)。 |
clear () |
|
copy () |
|
fromkeys (iterable[, value]) |
使用来自 iterable 的键创建一个新字典,并将值设置为 value。 |
get (key[, default]) |
如果字典中存在键 key,则返回其对应的值,否则返回默认值。 |
items () |
|
keys () |
|
pop (key[, default]) |
如果未找到 key,则如果提供了 default,则返回 default,否则引发 KeyError 异常。 |
popitem (/) |
移除并返回一个(key, value)对,作为一个二元组。 |
setdefault (key[, default]) |
如果 key 不在字典中,则将 key 插入,并将其值设置为 default。 |
update ([E, ]**F) |
如果 E 存在并且具有.keys()方法,则执行:对于 k 在 E 中:D[k] = E[k] 如果 E 存在但没有.keys()方法,则执行:对于 k, v 在 E 中:D[k] = v 在任一情况下,随后执行:对于 k 在 F 中:D[k] = F[k] |
values () |
scipy.optimize.OptimizeWarning
exception scipy.optimize.OptimizeWarning
with_traceback()
Exception.with_traceback(tb)
– 设置 self.traceback 为 tb 并返回 self。
scipy.optimize.minimize_scalar
scipy.optimize.minimize_scalar(fun, bracket=None, bounds=None, args=(), method=None, tol=None, options=None)
标量函数的局部最小化。
参数:
fun可调用对象
目标函数。标量函数,必须返回一个标量。
bracket序列,可选
对于方法‘brent’和‘golden’,bracket
定义了锁定间隔并且是必需的。可以是三元组(xa, xb, xc)
,满足xa < xb < xc
和func(xb) < func(xa) and func(xb) < func(xc)
,或者是一对(xa, xb)
,用于进行下山锁定搜索的初始点(参见scipy.optimize.bracket
)。最小化器res.x
不一定满足xa <= res.x <= xb
。
bounds序列,可选
对于方法‘bounded’,bounds是必须的,必须有两个有限项与优化边界相对应。
args元组,可选
传递给目标函数的额外参数。
method字符串或可调用对象,可选
求解器类型。应为以下之一:
- Brent
- Bounded
- Golden
- 自定义 - 可调用对象(从版本 0.14.0 开始添加),请参见下文
如果提供了边界,则默认为“Bounded”,否则为“Brent”。有关每个求解器的详细信息,请参见‘Notes’部分。
tol浮点数,可选
终止容差。要进行详细控制,请使用特定于解算器的选项。
options字典,可选
求解器选项的字典。
maxiterint
最大迭代次数。
dispbool
设置为 True 以打印收敛消息。
请参阅show_options
以获取特定于解算器的选项。
返回:
resOptimizeResult
表示优化结果的OptimizeResult
对象。重要属性包括:x
解决方案数组,success
指示优化器是否成功退出的布尔标志,message
描述终止原因。参见OptimizeResult
以获取其他属性的描述。
另请参见
用于标量多变量函数的最小化算法接口
解算器接受的额外选项
注意事项
本节描述了可以通过“method”参数选择的可用求解器。如果传递了bounds,默认方法是"Bounded"
Brent 方法;否则是无界的"Brent"
方法。
方法 Brent 使用 Brent 算法 [1] 寻找局部最小值。在可能的情况下,该算法使用反向抛物插值来加速黄金分割法的收敛速度。
方法 Golden 使用黄金分割搜索技术 [1]。它使用二分法的类似物来缩小括号内的区间。通常优先选择使用Brent方法。
方法 Bounded 可以执行有界最小化 [2] [3]。它使用 Brent 方法在区间 x1 < xopt < x2 中找到局部最小值。
注意,除非提供了有效的bracket
三元组,否则 Brent 和 Golden 方法不能保证成功。如果无法找到三点括号,请考虑使用scipy.optimize.minimize
。此外,所有方法仅用于局部最小化。当感兴趣的函数具有多个局部最小值时,请考虑全局优化。
自定义最小化器
当使用一些库的前端来进行minimize_scalar
时,传递自定义最小化方法可能很有用。您可以简单地将一个可调用对象作为method
参数传递。
可调用对象的调用形式为method(fun, args, **kwargs, **options)
,其中kwargs
对应于传递给minimize
的其他参数(如bracket
、tol等),除了options字典,其内容也会一对一地作为method参数传递。该方法应返回一个OptimizeResult
对象。
提供的method可调用对象必须能够接受(并可能忽略)任意参数;由于minimize
接受的参数集在将来版本中可能会扩展,这些参数也将一一传递给方法。您可以在 scipy.optimize 教程中找到一个例子。
版本 0.11.0 中的新功能。
参考文献
[1] (1,2)
Press, W., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes in C. Cambridge University Press.
[2]
Forsythe, G.E., M. A. Malcolm, and C. B. Moler. “Mathematical Computations 的计算机方法。” Prentice-Hall Series in Automatic Computation 259 (1977).
[3]
Brent, Richard P. Algorithms for Minimization Without Derivatives. Courier Corporation, 2013.
示例
考虑最小化以下函数的问题。
>>> def f(x):
... return (x - 2) * x * (x + 2)**2
使用Brent方法,我们找到了局部最小值如下:
>>> from scipy.optimize import minimize_scalar
>>> res = minimize_scalar(f)
>>> res.fun
-9.9149495908
最小化器是:
>>> res.x
1.28077640403
使用Bounded方法,我们找到了具有指定边界的局部最小值如下:
>>> res = minimize_scalar(f, bounds=(-3, -1), method='bounded')
>>> res.fun # minimum
3.28365179850e-13
>>> res.x # minimizer
-2.0000002026
scipy.optimize.minimize
scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)
最小化一个或多个变量的标量函数。
参数:
fun 可调用对象
要最小化的目标函数。
fun(x, *args) -> float
其中 x
是形状为 (n,) 的一维数组,args
是一个元组,包含完全指定函数所需的固定参数。
x0 数组,形状为 (n,)
初始猜测。大小为 (n,) 的实数元素数组,其中 n
是独立变量的数量。
args 元组,可选
传递给目标函数及其导数(fun、jac 和 hess 函数)的额外参数。
method 字符串或可调用对象,可选
求解器的类型。应为以下之一:
- ‘Nelder-Mead’ (详见此处)
- ‘Powell’ (详见此处)
- ‘CG’ (详见此处)
- ‘BFGS’ (详见此处)
- ‘Newton-CG’ (详见此处)
- ‘L-BFGS-B’ (详见此处)
- ‘TNC’ (详见此处)
- ‘COBYLA’ (详见此处)
- ‘SLSQP’ (详见此处)
- ‘trust-constr’(详见此处)
- ‘dogleg’ (详见此处)
- ‘trust-ncg’ (详见此处)
- ‘trust-exact’ (详见此处)
- ‘trust-krylov’ (详见此处)
- custom - 一个可调用对象,请参阅下文进行描述。
如果未提供,则根据问题是否有约束或边界选择 BFGS
、L-BFGS-B
、SLSQP
中的一种。
jac, optional
计算梯度向量的方法。仅适用于 CG、BFGS、Newton-CG、L-BFGS-B、TNC、SLSQP、dogleg、trust-ncg、trust-krylov、trust-exact 和 trust-constr。如果是可调用对象,则应为返回梯度向量的函数:
jac(x, *args) -> array_like, shape (n,)
其中 x
是形状为 (n,) 的数组,args
是具有固定参数的元组。如果 jac 是布尔值且为 True,则假定 fun 返回包含目标函数和梯度的元组 (f, g)
。方法 ‘Newton-CG’、‘trust-ncg’、‘dogleg’、‘trust-exact’ 和 ‘trust-krylov’ 要求提供一个可调用对象,或者 fun 返回目标函数和梯度。如果为 None 或 False,则使用绝对步长进行二点有限差分估计梯度。或者,关键字 {‘2-point’、‘3-point’、‘cs’} 可用于选择用于数值梯度估计的有限差分方案,并使用相对步长。这些有限差分方案遵循任何指定的 bounds。
hess,可选
计算 Hessian 矩阵的方法。仅适用于 Newton-CG、dogleg、trust-ncg、trust-krylov、trust-exact 和 trust-constr。如果是可调用对象,则应返回 Hessian 矩阵:
hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)
其中 x
是形状为 (n,) 的 ndarray,args
是具有固定参数的元组。关键字 {‘2-point’、‘3-point’、‘cs’} 也可用于选择用于数值估计 Hessian 的有限差分方案。或者,实现 HessianUpdateStrategy
接口的对象可用于近似 Hessian。实现此接口的可用拟牛顿方法包括:
并非每种方法都有所有选项;可参考注释中的可用性。
hessp可调用对象,可选
目标函数的 Hessian 矩阵乘以任意向量 p。仅适用于 Newton-CG、trust-ncg、trust-krylov、trust-constr。hessp 或 hess 之一需要提供。如果提供了 hess,则将忽略 hessp。hessp 必须计算 Hessian 乘以任意向量:
hessp(x, p, *args) -> ndarray 形状 (n,)
其中 x
是形状为 (n,) 的 ndarray,p
是维度为 (n,) 的任意向量,args
是具有固定参数的元组。
bounds序列或 Bounds
,可选
变量的界限用于 Nelder-Mead、L-BFGS-B、TNC、SLSQP、Powell、trust-constr 和 COBYLA 方法。有两种指定界限的方式:
Bounds
类的实例。- 对于 x 中的每个元素,
(min, max)
对用于指定界限。使用 None 表示无界限。
约束 或 约束列表 {约束,字典},可选
约束定义。仅适用于 COBYLA、SLSQP 和 trust-constr。
‘trust-constr’ 的约束被定义为一个单一对象或指定优化问题约束的对象列表。可用的约束类型包括:
LinearConstraint
NonlinearConstraint
对于 COBYLA、SLSQP,约束被定义为一个包含字段的字典列表:
typestr
约束类型:‘eq’ 表示等式约束,‘ineq’ 表示不等式约束。
funcallable
定义约束的函数。
jaccallable,可选
fun 的雅可比矩阵(仅适用于 SLSQP)。
argssequence,可选
传递给函数和雅可比矩阵的额外参数。
等式约束意味着约束函数的结果应为零,而不等式则意味着其应为非负。请注意,COBYLA 只支持不等式约束。
tolfloat,可选
终止容差。当指定 tol 时,所选的最小化算法设置一些相关的特定解算器容差为 tol。要进行详细的控制,请使用特定于解算器的选项。
optionsdict,可选
解算器选项的字典。除了 TNC 方法外,所有方法都接受以下通用选项:
maxiterint
执行的最大迭代次数。根据方法,每次迭代可能会使用多个函数评估。
对于 TNC 方法,请使用 maxfun 而不是 maxiter。
dispbool
设置为 True 以打印收敛消息。
对于特定方法的选项,请参见 show_options
。
callbackcallable,可选
每次迭代后调用的可调用函数。
除了 TNC、SLSQP 和 COBYLA 方法之外的所有方法都支持具有以下签名的可调用函数:
callback(intermediate_result: OptimizeResult)
其中 intermediate_result
是一个关键字参数,包含一个 OptimizeResult
,具有参数向量和目标函数当前值的属性。请注意,回调函数必须命名为 intermediate_result
,以便传递一个 OptimizeResult
。如果回调函数引发 StopIteration
,这些方法也将终止。
除了 trust-constr 方法之外的所有方法都支持以下形式的签名:
callback(xk)
其中 xk
是当前的参数向量。
使用内省来确定要调用上述哪种签名。
返回:
resOptimizeResult
优化结果表示为 OptimizeResult
对象。重要属性包括:x
解数组,success
表示优化器是否成功退出的布尔标志,message
描述终止原因。请参阅 OptimizeResult
了解其他属性的描述。
参见
minimize_scalar
标量函数最小化接口。
标量单变量函数最小化算法接口
show_options
函数显示选项。
求解器接受的额外选项
注意事项
本节描述可以通过 'method' 参数选择的可用求解器。默认方法为 BFGS。
无约束最小化
方法 CG 使用 Polak 和 Ribiere 的非线性共轭梯度算法,是 Fletcher-Reeves 方法的变种,详见 [5] pp.120-122. 仅使用一阶导数。
方法 BFGS 使用 Broyden、Fletcher、Goldfarb 和 Shanno(BFGS)拟牛顿法 [5] pp. 136. 仅使用一阶导数。即使在非平滑优化中,BFGS 也表现良好。此方法还返回存储在 OptimizeResult 对象的 hess_inv 中的海森矩阵逆的近似值。
方法 Newton-CG 使用 Newton-CG 算法 [5] pp. 168(也称为截断牛顿法)。它使用共轭梯度方法计算搜索方向。参见 TNC 方法,该方法类似,但适用于有边界约束的最小化问题。适用于大规模问题。
方法 dogleg 使用狗腿信任域算法 [5] 进行无约束最小化。该算法需要梯度和海森矩阵;此外,海森矩阵要求正定。
方法 trust-ncg 使用牛顿共轭梯度信任域算法 [5] 进行无约束最小化。该算法需要梯度和海森矩阵或计算给定向量与海森矩阵乘积的函数。适用于大规模问题。
方法 trust-krylov 使用 Newton GLTR 信赖域算法[14],[15]进行无约束最小化。此算法要求梯度和 Hessian 矩阵或计算给定向量与 Hessian 矩阵乘积的函数。适用于大规模问题。在不定问题上,通常比trust-ncg方法需要更少的迭代,推荐用于中等和大规模问题。
方法 trust-exact 是一种信赖域方法,用于无约束最小化,几乎完全解决二次子问题[13]。此算法要求梯度和 Hessian 矩阵(不要求为正定)。在许多情况下,这种方法收敛迭代较少,是小型和中型问题中最推荐的方法。
边界约束最小化
方法 Nelder-Mead 使用 Simplex 算法[1],[2]。此算法在许多应用中表现稳健。但是,如果可以信任数值导数的计算,其他利用一阶和/或二阶导数信息的算法可能更适合于其在一般情况下的更好性能。
方法 L-BFGS-B 使用 L-BFGS-B 算法[6],[7]进行边界约束最小化。
方法 Powell 是 Powell 方法的改进[3],[4],它是一种共轭方向方法。它沿着每个方向集合的每个向量(options和info中的direc字段)顺序进行一维最小化,每次主最小化循环迭代时更新。函数不需要可微,也不计算导数。如果未提供边界,则将使用无界线搜索。如果提供了边界,并且初始猜测在边界内,则最小化过程中的每个函数评估都将在边界内。如果提供了边界,初始猜测超出边界,并且direc具有完整秩(默认具有完整秩),则第一次迭代期间的某些函数评估可能超出边界,但第一次迭代后的每个函数评估都将在边界内。如果direc秩不完整,则某些参数可能不会被优化,并且不能保证解在边界内。
方法 TNC 使用截断牛顿算法[5], [8]来最小化带有变量界限的函数。此算法利用梯度信息;它也称为牛顿共轭梯度法。它与上述Newton-CG方法不同,因为它封装了 C 实现,并允许每个变量都有上限和下限。
约束最小化
方法 COBYLA 使用约束优化 BY 线性近似(COBYLA)方法[9], [10], [11]。该算法基于目标函数和每个约束的线性近似。该方法封装了该算法的 FORTRAN 实现。约束函数‘fun’可以返回单个数字或数字数组或列表。
方法 SLSQP 使用顺序最小二乘编程来最小化多变量函数,可以有各种边界、等式和不等式约束的组合。该方法封装了最初由 Dieter Kraft 实现的 SLSQP 优化子程序[12]。请注意,包装器通过将边界中的无限值转换为大浮点值来处理边界中的无限值。
方法 trust-constr 是一种用于约束优化的信赖域算法。它根据问题定义切换两种实现方式。这是 SciPy 中最通用的约束最小化算法,特别适用于大规模问题。对于等式约束问题,它是 Byrd-Omojokun 信赖域 SQP 方法的实现,详见[17]和[5],第 549 页。当还有不等式约束时,它切换到信赖域内点法,详见[16]。这种内点算法通过引入松弛变量并解决一系列逐渐减小的约束问题,以逐步减小的障碍参数。先前描述的等式约束 SQP 方法用于解决越来越精确的子问题,因为迭代逐渐接近解决方案。
有限差分期权
对于方法 trust-constr ,可以使用三种有限差分方案来近似梯度和海森矩阵:{‘2-point’, ‘3-point’, ‘cs’}。方案 ‘cs’ 可能是最精确的,但它要求函数能够正确处理复杂输入并在复平面上可微分。方案 ‘3-point’ 比 ‘2-point’ 更精确,但需要两倍的操作。如果通过有限差分法估计梯度,则必须使用一种拟牛顿策略来估计海森矩阵。
关于 hess 关键字的方法特定选项
method/Hess | None | callable | ‘2-point/’3-point’/’cs’ | HUS |
---|---|---|---|---|
Newton-CG | x | (n, n) LO | x | x |
dogleg | (n, n) | |||
trust-ncg | (n, n) | x | x | |
trust-krylov | (n, n) | x | x | |
trust-exact | (n, n) | |||
trust-constr | x | (n, n) LO sp | x | x |
其中 LO=LinearOperator,sp=Sparse matrix,HUS=HessianUpdateStrategy
自定义最小化器
可能在使用此方法的前端(例如scipy.optimize.basinhopping
)或其他库时,传递自定义的最小化方法可能很有用。您可以简单地将可调用对象作为method
参数传递。
可调用对象被调用为 method(fun, x0, args, **kwargs, **options)
,其中 kwargs
对应于传递给minimize
的任何其他参数(如 callback, hess 等),除了 options 字典,其内容也会逐对传递为 method 参数。此外,如果 jac 被传递为布尔类型,则 jac 和 fun 将被篡改,使 fun 仅返回函数值,而 jac 被转换为返回雅可比矩阵的函数。该方法应返回一个OptimizeResult
对象。
提供的 method 可调用对象必须能够接受(并可能忽略)任意参数;由于minimize
接受的参数集可能会在未来版本中扩展,这些参数将被传递给该方法。您可以在 scipy.optimize 教程中找到一个示例。
参考文献
[1]
Nelder, J A 和 R Mead. 1965. 函数最小化的单纯形法。《计算机期刊》 7: 308-13.
[2]
Wright M H. 1996. 直接搜索方法:曾经被蔑视,现在倍受尊重,收录于《数值分析 1995:1995 年邓迪双年会数值分析会议论文集》(主编 D F Griffiths 和 G A Watson)。Addison Wesley Longman, Harlow, UK. 191-208.
[3]
Powell, M J D. 1964. 一种在不计算导数的情况下找到多变量函数最小值的高效方法。《计算机期刊》 7: 155-162.
[4]
Press W, S A Teukolsky, W T Vetterling 和 B P Flannery. Numerical Recipes(任何版本),剑桥大学出版社。
[5] (1,2,3,4,5,6,7,8)
Nocedal, J 和 S J Wright. 2006.数值优化。Springer New York。
[6]
Byrd, R H 和 P Lu 和 J. Nocedal. 1995.用于有界约束优化的有限内存算法。SIAM Journal on Scientific and Statistical Computing 16(5):1190-1208。
[7]
Zhu, C 和 R H Byrd 和 J Nocedal. 1997. L-BFGS-B:算法 778:L-BFGS-B,FORTRAN 大规模有界约束优化的例程。ACM Transactions on Mathematical Software 23(4):550-560。
[8]
Nash, S G. 通过 Lanczos 方法的牛顿型最小化。1984. SIAM Journal of Numerical Analysis 21:770-778。
[9]
Powell, M J D. 一种直接搜索优化方法,通过线性插值模拟目标和约束函数。1994.优化和数值分析进展,主编 S. Gomez 和 J-P Hennart,Kluwer Academic(Dordrecht),51-67。
[10]
Powell M J D. 用于优化计算的直接搜索算法。1998. Acta Numerica 7:287-336。
[11]
Powell M J D. 无导数优化算法概览。2007.剑桥大学技术报告 DAMTP 2007/NA03
[12]
Kraft, D. 用于顺序二次规划的软件包。1988. Tech. Rep. DFVLR-FB 88-28,DLR German Aerospace Center – Institute for Flight Mechanics,Koln,Germany。
[13]
Conn, A. R., Gould, N. I.,和 Toint, P. L. 信任区域方法。2000. Siam. pp. 169-200.
[14]
F. Lenders, C. Kirches, A. Potschka: “trlib:用于迭代解决信任区域问题的无向量实现”,arXiv:1611.04718
[15]
N. Gould, S. Lucidi, M. Roma, P. Toint: “使用 Lanczos 方法解决信任区域子问题”,SIAM J. Optim., 9(2), 504–525, (1999).
[16]
Byrd, Richard H., Mary E. Hribar 和 Jorge Nocedal. 1999.大规模非线性规划的内点算法。SIAM Journal on Optimization 9.4:877-900。
[17]
Lalee, Marucha,Jorge Nocedal 和 Todd Plantega. 1998.关于大规模等式约束优化算法的实现。SIAM Journal on Optimization 8.3:682-706。
例子
让我们考虑最小化 Rosenbrock 函数的问题。该函数(及其相应的导数)在rosen
(分别在rosen_der
,rosen_hess
中实现)中。scipy.optimize
。
>>> from scipy.optimize import minimize, rosen, rosen_der
Nelder-Mead方法的一个简单应用是:
>>> x0 = [1.3, 0.7, 0.8, 1.9, 1.2]
>>> res = minimize(rosen, x0, method='Nelder-Mead', tol=1e-6)
>>> res.x
array([ 1., 1., 1., 1., 1.])
现在使用BFGS算法,使用第一阶导数和一些选项:
>>> res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
... options={'gtol': 1e-6, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 26
Function evaluations: 31
Gradient evaluations: 31
>>> res.x
array([ 1., 1., 1., 1., 1.])
>>> print(res.message)
Optimization terminated successfully.
>>> res.hess_inv
array([
[ 0.00749589, 0.01255155, 0.02396251, 0.04750988, 0.09495377], # may vary
[ 0.01255155, 0.02510441, 0.04794055, 0.09502834, 0.18996269],
[ 0.02396251, 0.04794055, 0.09631614, 0.19092151, 0.38165151],
[ 0.04750988, 0.09502834, 0.19092151, 0.38341252, 0.7664427 ],
[ 0.09495377, 0.18996269, 0.38165151, 0.7664427, 1.53713523]
])
接下来,考虑一个带有多个约束条件的最小化问题(即来自[5]的示例 16.4)。目标函数是:
>>> fun = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2
有三个定义为约束条件:
>>> cons = ({'type': 'ineq', 'fun': lambda x: x[0] - 2 * x[1] + 2},
... {'type': 'ineq', 'fun': lambda x: -x[0] - 2 * x[1] + 6},
... {'type': 'ineq', 'fun': lambda x: -x[0] + 2 * x[1] + 2})
变量必须为正数,因此以下是界限:
>>> bnds = ((0, None), (0, None))
优化问题使用 SLSQP 方法求解如下:
>>> res = minimize(fun, (2, 0), method='SLSQP', bounds=bnds,
... constraints=cons)
它应该收敛到理论解(1.4, 1.7)。
scipy.optimize.NonlinearConstraint
class scipy.optimize.NonlinearConstraint(fun, lb, ub, jac='2-point', hess=<scipy.optimize._hessian_update_strategy.BFGS object>, keep_feasible=False, finite_diff_rel_step=None, finite_diff_jac_sparsity=None)
变量的非线性约束。
约束具有一般的不等式形式:
lb <= fun(x) <= ub
这里独立变量向量 x 作为形状为 (n,) 的 ndarray 传递,并且fun
返回一个有 m 个分量的向量。
可以使用相等的边界来表示等式约束,或者使用无限边界来表示单边约束。
参数:
fun:可调用函数
定义约束的函数。签名是 fun(x) -> array_like, shape (m,)
。
lb, ub:array_like
约束的下界和上界。每个数组必须具有形状 (m,) 或者是一个标量,后一种情况下约束对所有约束分量是相同的。使用 np.inf
与适当的符号指定单边约束。将 lb 和 ub 的分量设置为相等表示等式约束。注意,可以通过设置 lb 和 ub 的不同分量来混合不同类型的约束:区间、单边或等式。
jac:{可调用函数, ‘2-point’, ‘3-point’, ‘cs’},可选
计算雅可比矩阵的方法(一个 m×n 矩阵,其中元素 (i, j) 是 f[i] 对 x[j] 的偏导数)。关键字 {‘2-point’, ‘3-point’, ‘cs’} 选择数值估计的有限差分方案。一个可调用对象必须具有以下签名:jac(x) -> {ndarray, sparse matrix}, shape (m, n)
。默认为 ‘2-point’。
hess:{可调用函数, ‘2-point’, ‘3-point’, ‘cs’, HessianUpdateStrategy, None},可选
计算 Hessian 矩阵的方法。关键字 {‘2-point’, ‘3-point’, ‘cs’} 选择数值估计的有限差分方案。或者,实现了 HessianUpdateStrategy
接口的对象可以用来近似 Hessian。当前可用的实现是:
一个可调用对象必须返回 dot(fun, v)
的 Hessian 矩阵,并且必须具有以下签名:hess(x, v) -> {LinearOperator, sparse matrix, array_like}, shape (n, n)
。这里 v
是形状为 (m,) 的 ndarray,包含拉格朗日乘数。
keep_feasible:bool 类型的数组,可选
决定在迭代过程中是否保持约束分量的可行性。一个单一的值设置该属性对所有分量生效。默认值为 False。对于等式约束没有影响。
finite_diff_rel_step:None 或者 array_like,可选
有限差分近似的相对步长。默认为 None,根据有限差分方案自动选择合理值。
finite_diff_jac_sparsity: {None, array_like, sparse matrix}, optional
定义了有限差分估计雅可比矩阵的稀疏结构,其形状必须为 (m, n)。如果雅可比矩阵每行只有少量非零元素,在提供稀疏结构的情况下将大大加快计算速度。零条目意味着雅可比矩阵中对应元素恒为零。如果提供,则强制使用 'lsmr' 信赖区域求解器。如果为 None(默认值),则将使用稠密差分。
注意事项
可用于近似雅可比或海森矩阵的有限差分方案 {‘2-point’, ‘3-point’, ‘cs’}。然而,我们不允许同时用于两者的近似。因此,每当通过有限差分估计雅可比时,我们要求使用一种拟牛顿策略估计海森矩阵。
方案 'cs' 可能是最准确的,但要求函数能正确处理复杂输入,并在复平面上解析延拓。方案 '3-point' 比 '2-point' 更精确,但操作数量是其两倍。
示例
约束条件 x[0] < sin(x[1]) + 1.9
>>> from scipy.optimize import NonlinearConstraint
>>> import numpy as np
>>> con = lambda x: x[0] - np.sin(x[1])
>>> nlc = NonlinearConstraint(con, -np.inf, 1.9)
scipy.optimize.LinearConstraint
class scipy.optimize.LinearConstraint(A, lb=-inf, ub=inf, keep_feasible=False)
变量的线性约束。
约束具有一般不等式形式:
lb <= A.dot(x) <= ub
这里作为独立变量向量 x 以形状为(m,)的 ndarray 传递,矩阵 A 的形状为(m, n)。
可以使用相等的边界来表示等式约束或无穷边界来表示单侧约束。
参数:
A,形状为(m, n)
定义约束的矩阵。
lb, ub稠密的数组,可选
约束的下限和上限。每个数组必须具有形状(m,)或者是标量,在后一种情况下,约束的所有组件都将具有相同的边界。使用np.inf
和适当的符号来指定单侧约束。将lb和ub的组件设置相等以表示等式约束。请注意,可以通过根据需要设置lb和ub的不同组件来混合不同类型的约束:区间约束、单侧约束或等式约束。默认为lb = -np.inf
和ub = np.inf
(无限制)。
keep_feasible稠密的布尔数组,可选
是否在迭代过程中保持约束组件的可行性。单个值设置此属性以适用于所有组件。默认为 False。对于等式约束没有影响。
方法
residual (x) |
计算约束函数与限制之间的残差 |
---|
scipy.optimize.Bounds
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.Bounds.html#scipy.optimize.Bounds
class scipy.optimize.Bounds(lb=-inf, ub=inf, keep_feasible=False)
变量的边界约束。
约束条件具有一般的不等式形式:
lb <= x <= ub
可以使用相等的边界表示等式约束或无穷大的边界表示单侧约束。
参数:
lb, ub稠密的数组,可选
自变量的下限和上限。lb、ub和keep_feasible必须具有相同的形状或可广播。将lb和ub的组件设为相等以固定变量。使用np.inf
和适当的符号禁用所有或部分变量的边界。请注意,可以通过必要时设置lb和ub的不同组件来混合不同类型的约束:区间约束、单侧约束或等式约束。默认为lb = -np.inf
和ub = np.inf
(无边界)。
keep_feasible稠密的 bool 数组,可选
是否在迭代过程中保持约束组件的可行性。必须与lb和ub进行广播。默认为 False。对于等式约束无影响。
方法
residual (x) |
计算输入与边界之间的残差(松弛度) |
---|
scipy.optimize.BFGS
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.BFGS.html#scipy.optimize.BFGS
class scipy.optimize.BFGS(exception_strategy='skip_update', min_curvature=None, init_scale='auto')
Broyden-Fletcher-Goldfarb-Shanno (BFGS) Hessian 更新策略。
参数:
exception_strategy, 可选
定义在曲率条件违反时如何进行。将其设置为 ‘skip_update’ 以跳过更新。或者,将其设置为 ‘damp_update’ 以在实际的 BFGS 结果和未修改的矩阵之间插值。这两种异常策略在 [1],p.536-537 中有解释。
min_curvaturefloat
该数字乘以归一化因子,定义了允许不受异常策略影响的最小曲率 dot(delta_grad, delta_x)
。当 exception_strategy = 'skip_update'
时,默认为 1e-8,当 exception_strategy = 'damp_update'
时,默认为 0.2。
init_scale
矩阵在第一次迭代时的尺度。在第一次迭代中,Hessian 矩阵或其逆将初始化为 init_scale*np.eye(n)
,其中 n
是问题的维度。将其设置为 ‘auto’ 可以使用自动启发式方法选择初始尺度。该启发式方法在 [1],p.143 中描述。默认使用 ‘auto’。
注意事项
更新基于 [1],p.140 中的描述。
参考文献
[1] (1,2,3)
Nocedal, Jorge, and Stephen J. Wright. “数值优化” 第二版 (2006)。
方法
dot (p) |
计算内部矩阵与给定向量的乘积。 |
---|---|
get_matrix () |
返回当前内部矩阵。 |
initialize (n, approx_type) |
初始化内部矩阵。 |
update (delta_x, delta_grad) |
更新内部矩阵。 |
scipy.optimize.SR1
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.SR1.html#scipy.optimize.SR1
class scipy.optimize.SR1(min_denominator=1e-08, init_scale='auto')
对称秩 1 Hessian 更新策略。
参数:
min_denominator浮点数
此数字通过归一化因子缩放,定义了更新中允许的最小分母大小。当条件违反时,我们会跳过更新。默认使用1e-8
。
init_scale,可选
在第一次迭代中,Hessian 矩阵或其逆将用init_scale*np.eye(n)
初始化,其中n
是问题的维度。将其设置为'auto',以便使用自动启发式方法选择初始规模。该启发式方法在[1],p.143 中描述。默认情况下使用'auto'。
注意
更新基于描述[1],p.144-146。
参考文献
[1] (1,2)
Nocedal, Jorge, and Stephen J. Wright. “Numerical optimization” Second Edition (2006).
方法
dot (p) |
计算内部矩阵与给定向量的乘积。 |
---|---|
get_matrix () |
返回当前内部矩阵。 |
initialize (n, approx_type) |
初始化内部矩阵。 |
update (delta_x, delta_grad) |
更新内部矩阵。 |
scipy.optimize.basinhopping
scipy.optimize.basinhopping(func, x0, niter=100, T=1.0, stepsize=0.5, minimizer_kwargs=None, take_step=None, accept_test=None, callback=None, interval=50, disp=False, niter_success=None, seed=None, *, target_accept_rate=0.5, stepwise_factor=0.9)
使用盆地跳跃算法寻找函数的全局最小值。
盆地跳跃是一种两阶段方法,结合了全局步进算法和每步的局部最小化。设计用于模仿原子簇能量最小化的自然过程,适用于具有“漏斗形但崎岖”的能量景观的类似问题 [5]。
由于步骤采取、步骤接受和最小化方法都是可定制的,因此该函数也可以用于实现其他两阶段方法。
参数:
funccallable f(x, *args)
要优化的函数。args
可以作为字典 minimizer_kwargs 的可选项传递。
x0array_like
初始猜测。
niter整数,可选
盆地跳跃迭代次数。将有 niter + 1
次局部最小化运行。
T浮点数,可选
接受或拒绝标准的“温度”参数。较高的“温度”意味着将接受更大的函数值跳跃。为了获得最佳结果,T 应与局部最小值之间的分离(在函数值上)相当。
stepsize浮点数,可选
用于随机位移的最大步长。
minimizer_kwargsdict,可选
要传递给本地最小化器 scipy.optimize.minimize
的额外关键字参数。一些重要的选项可能包括:
methodstr
最小化方法(例如
"L-BFGS-B"
)argstuple
传递给目标函数 (func) 及其导数(Jacobian、Hessian)的额外参数。
take_stepcallable take_step(x)
,可选
用此例程替换默认的步进例程。默认的步进例程是坐标的随机位移,但其他步进算法可能对某些系统更好。take_step 可以选择具有属性 take_step.stepsize
。如果存在此属性,则 basinhopping
将调整 take_step.stepsize
以尝试优化全局最小搜索。
accept_testcallable,accept_test(f_new=f_new, x_new=x_new, f_old=fold, x_old=x_old)
,可选
定义一个测试,用于判断是否接受该步骤。这将与基于“温度”T的 Metropolis 测试一起使用。可以接受的返回值为 True、False 或"force accept"
。如果任何测试返回 False,则拒绝该步骤。如果是后者,则这将覆盖任何其他测试以接受该步骤。例如,可以强制性地从basinhopping
被困住的局部最小值中逃脱。
callback可调用对象,callback(x, f, accept)
,可选
一个回调函数,将为找到的所有最小值调用。x和f是试探最小值的坐标和函数值,accept表示是否接受该最小值。例如,可以使用此功能保存找到的最低的 N 个最小值。此外,callback可以用于通过可选地返回 True 来指定用户定义的停止标准,以停止basinhopping
程序运行。
interval整数,可选
用于定期更新stepsize的间隔
disp布尔值,可选
设置为 True 以打印状态消息
niter_success整数,可选
如果全局最小候选在此迭代次数内保持不变,则停止运行。
seed{None, 整数, numpy.random.Generator
, numpy.random.RandomState
},可选
如果seed为 None(或np.random),则使用numpy.random.RandomState
单例。如果seed是整数,则使用一个新的RandomState
实例,并使用seed进行种子化。如果seed已经是Generator
或RandomState
实例,则直接使用该实例。指定seed以进行可重复的最小化。使用此种子生成的随机数仅影响默认的 Metropolis accept_test和默认的take_step。如果您提供自己的take_step和accept_test,并且这些函数使用随机数生成,则这些函数负责其随机数生成器的状态。
target_accept_rate浮点数,可选
用于调整stepsize的目标接受率。如果当前接受率大于目标,则增加stepsize。否则,减少stepsize。范围为(0, 1)。默认值为 0.5。
自版本 1.8.0 新增。
stepwise_factor浮点数,可选
stepsize在每次更新时乘以或除以此步进因子。范围为(0, 1)。默认值为 0.9。
自版本 1.8.0 新增。
返回:
resOptimizeResult
优化结果表示为OptimizeResult
对象。重要属性包括:x
解数组,fun
解处函数值,以及 message
描述终止原因。所选最小化器在最低最小值处返回的OptimizeResult
对象也包含在此对象中,并可通过lowest_optimization_result
属性访问。参见OptimizeResult
获取其他属性的描述。
另见
minimize
局部最小化函数对每个 basinhopping 步骤调用。minimizer_kwargs 被传递给此例程。
注意事项
Basin-hopping 是一种随机算法,旨在找到一个或多个变量的光滑标量函数的全局最小值[1] [2] [3] [4]。该算法在目前的形式下由 David Wales 和 Jonathan Doye 描述[2] www-wales.ch.cam.ac.uk/
。
算法是迭代的,每个周期由以下特征组成
-
坐标的随机扰动
-
局部最小化
-
基于最小化函数值接受或拒绝新坐标
此处使用的接受测试是标准 Monte Carlo 算法的 Metropolis 标准,尽管还有许多其他可能性[3]。
已证明该全局最小化方法对物理和化学中的各种问题非常高效。当函数具有由大障碍物分隔的多个极小值时特别有用。请参见Cambridge Cluster Database以获取主要使用 basin-hopping 优化的分子系统数据库。该数据库包括超过 300 个自由度的最小化问题。
有关 basin-hopping 的 Fortran 实现请参见自由软件程序GMIN。该实现包括上述过程的许多变体,包括更高级的步骤算法和替代接受标准。
对于随机全局优化,无法确定是否实际上找到了真正的全局最小值。作为一致性检查,可以从多个不同的随机起始点运行算法,以确保每个示例中找到的最低最小值已收敛到全局最小值。因此,默认情况下,basinhopping
将仅运行 niter 次迭代,并返回找到的最低最小值。用户需要确保这实际上是全局最小值。
选择 stepsize:这是 basinhopping
中的关键参数,取决于正在解决的问题。步长在每个维度中均匀选择,从 x0-stepsize 到 x0+stepsize 的区域内。理想情况下,它应与被优化函数的局部极小值之间(在参数值上的)典型分离可比较。basinhopping
将默认调整 stepsize 以找到最优值,但这可能需要多次迭代。如果设置一个合理的初始值给 stepsize
,则可以更快地获得结果。
选择 T:参数 T 是 Metropolis 准则中使用的“温度”。如果 func(xnew) < func(xold)
,则 Basin-hopping 步骤始终被接受。否则,它们将以以下概率被接受:
exp( -(func(xnew) - func(xold)) / T )
因此,为了获得最佳结果,T 应与局部极小值之间(在函数值上的)典型差异可比较。(“墙”高度对局部极小值无关紧要。)
如果 T 为 0,则算法变为单调 Basin-Hopping,其中所有增加能量的步骤都被拒绝。
0.12.0 版本中的新内容。
参考文献
[1]
Wales, David J. 2003, 能量景观,剑桥大学出版社,剑桥,英国。
[2] (1,2)
Wales, D J 和 Doye J P K, Lennard-Jones 簇的基态结构的全局优化:通过 Basin-Hopping 和包含多达 110 个原子的结构。《物理化学学报》,1997 年,101,5111。
[3] (1,2)
Li, Z. 和 Scheraga, H. A., 蛋白质折叠中的多极小问题的蒙特卡洛-最小化方法,《美国国家科学院院刊》,1987 年,84,6611。
[4]
Wales, D. J. 和 Scheraga, H. A., 簇、晶体和生物分子的全局优化,《科学》,1999 年,285,1368。
[5]
Olson, B., Hashmi, I., Molloy, K., 和 Shehu1, A., Basin Hopping 作为生物大分子特征化的一般和多功能优化框架,《人工智能进展》,2012 年卷(2012),文章 ID 674832,DOI:10.1155/2012/674832
例子
下面的例子是一个一维最小化问题,在抛物线上叠加了许多局部极小值。
>>> import numpy as np
>>> from scipy.optimize import basinhopping
>>> func = lambda x: np.cos(14.5 * x - 0.3) + (x + 0.2) * x
>>> x0 = [1.]
盆地跳跃内部使用局部最小化算法。我们将使用参数minimizer_kwargs告诉盆地跳跃使用哪种算法以及如何设置该最小化器。此参数将传递给scipy.optimize.minimize
。
>>> minimizer_kwargs = {"method": "BFGS"}
>>> ret = basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs,
... niter=200)
>>> print("global minimum: x = %.4f, f(x) = %.4f" % (ret.x, ret.fun))
global minimum: x = -0.1951, f(x) = -1.0009
接下来考虑一个二维最小化问题。此外,这次我们将使用梯度信息来显著加速搜索。
>>> def func2d(x):
... f = np.cos(14.5 * x[0] - 0.3) + (x[1] + 0.2) * x[1] + (x[0] +
... 0.2) * x[0]
... df = np.zeros(2)
... df[0] = -14.5 * np.sin(14.5 * x[0] - 0.3) + 2. * x[0] + 0.2
... df[1] = 2. * x[1] + 0.2
... return f, df
我们还将使用不同的局部最小化算法。此外,我们必须告诉最小化器,我们的函数同时返回能量和梯度(雅可比矩阵)。
>>> minimizer_kwargs = {"method":"L-BFGS-B", "jac":True}
>>> x0 = [1.0, 1.0]
>>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
... niter=200)
>>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0],
... ret.x[1],
... ret.fun))
global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109
下面是一个使用自定义步进例程的示例。想象一下,你希望第一个坐标采取比其他坐标更大的步骤。可以这样实现:
>>> class MyTakeStep:
... def __init__(self, stepsize=0.5):
... self.stepsize = stepsize
... self.rng = np.random.default_rng()
... def __call__(self, x):
... s = self.stepsize
... x[0] += self.rng.uniform(-2.*s, 2.*s)
... x[1:] += self.rng.uniform(-s, s, x[1:].shape)
... return x
由于MyTakeStep.stepsize
存在,盆地跳跃将调整stepsize的大小以优化搜索。我们将使用与之前相同的二维函数。
>>> mytakestep = MyTakeStep()
>>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
... niter=200, take_step=mytakestep)
>>> print("global minimum: x = [%.4f, %.4f], f(x) = %.4f" % (ret.x[0],
... ret.x[1],
... ret.fun))
global minimum: x = [-0.1951, -0.1000], f(x) = -1.0109
现在,让我们使用一个自定义回调函数的示例,该函数打印出每个找到的最小值的值。
>>> def print_fun(x, f, accepted):
... print("at minimum %.4f accepted %d" % (f, int(accepted)))
这次我们只运行 10 次盆地跳步。
>>> rng = np.random.default_rng()
>>> ret = basinhopping(func2d, x0, minimizer_kwargs=minimizer_kwargs,
... niter=10, callback=print_fun, seed=rng)
at minimum 0.4159 accepted 1
at minimum -0.4317 accepted 1
at minimum -1.0109 accepted 1
at minimum -0.9073 accepted 1
at minimum -0.4317 accepted 0
at minimum -0.1021 accepted 1
at minimum -0.7425 accepted 1
at minimum -0.9073 accepted 1
at minimum -0.4317 accepted 0
at minimum -0.7425 accepted 1
at minimum -0.9073 accepted 1
在第 8 次迭代中已经找到的最小值为-1.0109,实际上是全局最小值。
scipy.optimize.brute
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.brute.html#scipy.optimize.brute
scipy.optimize.brute(func, ranges, args=(), Ns=20, full_output=0, finish=<function fmin>, disp=False, workers=1)
通过蛮力法在给定范围内最小化一个函数。
使用“蛮力”方法,即在多维点网格的每个点计算函数的值,以找到函数的全局最小值。
函数在调用时以第一个调用函数的数据类型在范围内进行评估,由 vectorize
NumPy 函数强制执行。当 full_output=True
时,函数评估的值和类型受 finish 参数的影响(详见 Notes)。
蛮力法是低效的,因为网格点的数量呈指数增长 - 要评估的网格点数量为 Ns ** len(x)
。因此,即使是粗略的网格间距,中等规模的问题也可能需要很长时间运行,或者会遇到内存限制。
参数:
funccallable
要最小化的目标函数。必须是形式为 f(x, *args)
的函数,其中 x
是一个一维数组的参数,而 args
是一个元组,包含完全指定函数所需的任何额外固定参数。
rangestuple
ranges 元组的每个组件必须是“切片对象”或形如 (low, high)
的范围元组。程序使用这些来创建网格点,以便计算目标函数。详见 Note 2。
argstuple, optional
任何额外固定参数,以完全指定函数。
Nsint, optional
如果未另有说明,每个轴上的网格点数。详见 Note2。
full_outputbool, optional
如果为 True,则返回评估网格及其上的目标函数值。
finishcallable, optional
一个优化函数,它以蛮力最小化的结果作为初始猜测进行调用。finish 应将 func 和初始猜测作为位置参数,并将 args 作为关键字参数。它还可以作为关键字参数接受 full_output 和/或 disp。如果不使用“抛光”函数,则使用 None。详见 Notes 获取更多详情。
dispbool, optional
设置为 True 时,打印来自 finish 可调用的收敛消息。
workersint or map-like callable, optional
如果 workers 是一个整数,则将网格细分为 workers 部分,并并行评估(使用 multiprocessing.Pool
)。提供 -1 使用所有可用的核心进程。或者提供一个类似映射的可调用对象,例如 multiprocessing.Pool.map 用于并行评估网格。此评估是作为 workers(func, iterable)
进行的。要求 func 可被 pickle。
自版本 1.3.0 起新增。
返回:
x0ndarray
包含目标函数取得其最小值的点的坐标的一维数组。(参见注 1,了解返回的是哪个点。)
fval浮点数
x0点处的函数值。(当full_output为 True 时返回。)
grid元组
评估网格的表示。它与x0的长度相同。(当full_output为 True 时返回。)
Jout数组
在评估网格的每个点处的函数值,即Jout = func(*grid)
。(当full_output为 True 时返回。)
另请参见
basinhopping
,differential_evolution
笔记
注 1:程序找到了目标函数取得最低值的网格点。如果finish为 None,则返回该点。当全局最小值出现在网格边界内(或者非常接近)且网格足够精细时,该点将位于全局最小值的邻近区域。
然而,用户通常使用其他优化程序来“磨光”网格点的值,即在brute的最佳网格点附近寻找更精确(局部)的最小值。brute
函数的finish选项提供了一种方便的方法来实现此目的。使用的任何磨光程序必须将brute的输出作为其位置参数的初始猜测,并将brute的输入值作为关键字参数的args。否则将会引发错误。它还可以作为关键字参数接受full_output和/或disp。
brute
假设finish函数返回OptimizeResult
对象或形如(xmin, Jmin, ... , statuscode)
的元组,其中xmin
是参数的最小值,Jmin
是目标函数的最小值,“…”可能是其他返回的值(brute
不使用),而statuscode
是finish程序的状态码。
注意,当finish不为 None 时,返回的值是finish程序的结果,而不是网格点的结果。因此,虽然brute
限制其搜索在输入网格点上,但finish程序的结果通常不会与任何网格点重合,并且可能落在网格的边界之外。因此,如果仅需要在提供的网格点上找到最小值,请确保传入finish=None。
注 2:点的网格是一个numpy.mgrid
对象。对于brute
,ranges和Ns的输入具有以下效果。ranges元组的每个组件可以是一个切片对象或一个给定值范围的两元组,比如(0, 5)。如果组件是一个切片对象,brute
直接使用它。如果组件是一个两元组范围,brute
内部将其转换为一个切片对象,该对象从其低值到其高值插值出Ns个点,包括两端的值。
例子
我们演示了使用brute
来寻找一个由正定二次型和两个深“高斯形”坑的函数的全局最小值。具体地,定义目标函数f为另外三个函数的和,f = f1 + f2 + f3
。我们假设每个函数都有一个签名(z, *params)
,其中z = (x, y)
,而params
和函数如下所定义。
>>> import numpy as np
>>> params = (2, 3, 7, 8, 9, 10, 44, -1, 2, 26, 1, -2, 0.5)
>>> def f1(z, *params):
... x, y = z
... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
... return (a * x**2 + b * x * y + c * y**2 + d*x + e*y + f)
>>> def f2(z, *params):
... x, y = z
... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
... return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
>>> def f3(z, *params):
... x, y = z
... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
... return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
>>> def f(z, *params):
... return f1(z, *params) + f2(z, *params) + f3(z, *params)
因此,目标函数可能在由其组成的三个函数的最小值附近有局部极小值。为了使用fmin
来优化其格点结果,我们可以继续如下操作:
>>> rranges = (slice(-4, 4, 0.25), slice(-4, 4, 0.25))
>>> from scipy import optimize
>>> resbrute = optimize.brute(f, rranges, args=params, full_output=True,
... finish=optimize.fmin)
>>> resbrute[0] # global minimum
array([-1.05665192, 1.80834843])
>>> resbrute[1] # function value at global minimum
-3.4085818767
请注意,如果finish被设置为 None,我们将得到格点[-1.0 1.75],其中四舍五入的函数值为-2.892。
scipy.optimize.differential_evolution
scipy.optimize.differential_evolution(func, bounds, args=(), strategy='best1bin', maxiter=1000, popsize=15, tol=0.01, mutation=(0.5, 1), recombination=0.7, seed=None, callback=None, disp=False, polish=True, init='latinhypercube', atol=0, updating='immediate', workers=1, constraints=(), x0=None, *, integrality=None, vectorized=False)
找到多元函数的全局最小值。
差分进化方法[1]具有随机性质。它不使用梯度方法来找到最小值,可以搜索候选空间的大面积,但通常需要比传统的基于梯度的技术更多的函数评估次数。
算法是由斯托恩和普赖斯[2]提出的。
参数:
func可调用
要最小化的目标函数。必须以 f(x, *args)
的形式存在,其中 x 是形式为 1-D 数组的参数,args 是任何额外的固定参数元组,用于完全指定函数。参数数量 N 等于 len(x)
。
bounds序列或Bounds
变量的边界。有两种指定边界的方式:
Bounds
类的实例。(min, max)
对,用于定义 func 优化参数 x 的有限下限和上限。
使用总边界数来确定参数数量 N。如果有参数的边界相等,则自由参数的总数为 N - N_equal
。
args元组,可选
用于完全指定目标函数的任何额外固定参数。
strategy,可选
使用的差分进化策略。应为以下之一:
- ‘best1bin’
- ‘best1exp’
- ‘rand1bin’
- ‘rand1exp’
- ‘rand2bin’
- ‘rand2exp’
- ‘randtobest1bin’
- ‘randtobest1exp’
- ‘currenttobest1bin’
- ‘currenttobest1exp’
- ‘best2exp’
- ‘best2bin’
默认为‘best1bin’。可以实施的策略在‘Notes’中有概述。另外,差分进化策略可以通过提供一个构建试验向量的可调用对象来进行定制化。该可调用对象必须具有形式strategy(candidate: int, population: np.ndarray, rng=None)
,其中candidate
是一个整数,指定要进化的种群条目,population
是形状为(S, N)
的数组,包含所有种群成员(其中 S 是总种群大小),rng
是求解器内使用的随机数生成器。candidate
的范围为[0, S)
。strategy
必须返回一个形状为(N,)的试验向量。将此试验向量的适应度与population[candidate]
的适应度进行比较。
自版本 1.12.0 起更改:通过可调用对象定制进化策略。
maxiterint,可选
演化整个种群的最大代数。没有优化的最大函数评估次数为:(maxiter + 1) * popsize * (N - N_equal)
。
popsizeint, optional
用于设置总种群大小的乘数。种群包含popsize * (N - N_equal)
个个体。如果通过init关键字提供了初始种群,则此关键字将被覆盖。使用init='sobol'
时,种群大小计算为popsize * (N - N_equal)
后的下一个 2 的幂。
tolfloat, optional
收敛的相对容差,当np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))
时停止求解,其中atol和tol分别为绝对容差和相对容差。
mutationfloat or tuple(float, float), optional
变异常数。在文献中,这也被称为差分权重,用 F 表示。如果指定为浮点数,则应在[0, 2]范围内。如果指定为元组(min, max)
,则采用抖动。抖动会随机地在每代基础上改变变异常数。该代的变异常数取自U[min, max)
。抖动可以显著加快收敛速度。增加变异常数会增加搜索半径,但会减慢收敛速度。
recombinationfloat, optional
重组常数,应在[0, 1]范围内。在文献中,这也被称为交叉概率,用 CR 表示。增加此值允许更多的突变体进入下一代,但会有种群稳定性的风险。
seed{None, int, numpy.random.Generator
, numpy.random.RandomState
}, optional
如果seed为 None(或np.random),则使用numpy.random.RandomState
单例。如果seed为整数,则使用一个新的RandomState
实例,并用seed进行种子化。如果seed已经是Generator
或RandomState
实例,则使用该实例。为了重复最小化操作,请指定seed。
dispbool, optional
每次迭代时打印评估的func。
callbackcallable, optional
每次迭代后调用的可调用对象,具有以下签名:
callback(intermediate_result: OptimizeResult)
。
其中intermediate_result
是一个包含OptimizeResult
属性x
和fun
的关键字参数,表示到目前为止找到的最佳解和目标函数值。请注意,回调函数必须传递一个名为intermediate_result
的OptimizeResult
。
回调还支持类似的签名:
callback(x, convergence: float=val)
val
表示种群收敛的分数值。当val
大于1.0
时,函数停止。
使用内省来确定调用的签名之一。
全局最小化将在回调引发StopIteration
或返回True
时终止;仍会执行任何调整。
在版本 1.12.0 中更改:回调接受intermediate_result
关键字。
polishbool,可选
如果为 True(默认),则最后使用L-BFGS-B方法的scipy.optimize.minimize
来优化最佳种群成员,这可能会略微改善最小化。如果正在研究受约束问题,则改为使用trust-constr方法。对于具有许多约束的大问题,由于雅可比计算,优化可能需要很长时间。
initstr 或类数组,可选
指定执行的种群初始化类型。应为以下之一:
- ‘latinhypercube’
- ‘sobol’
- ‘halton’
- ‘random’
- 数组指定初始种群。数组应具有形状
(S, N)
,其中 S 为总种群大小,N 为参数数量。在使用之前,init将被剪辑到bounds内。
默认为‘latinhypercube’。拉丁超立方体采样试图最大化可用参数空间的覆盖。
‘sobol’和‘halton’是更优的选择,甚至更大程度地最大化参数空间。‘sobol’将强制初始种群大小计算为popsize * (N - N_equal)
后的下一个 2 的幂。‘halton’没有要求,但效率稍低。有关更多详情,请参见scipy.stats.qmc
。
‘random’随机初始化种群 - 这样做的缺点是可能会发生聚类,从而阻止参数空间的整体覆盖。例如,可以使用数组来指定种群,以在已知解存在的位置创建一组紧密的初始猜测,从而减少收敛时间。
atolfloat,可选
收敛的绝对容差,当np.std(pop) <= atol + tol * np.abs(np.mean(population_energies))
时,求解停止,其中atol和tol分别是绝对容差和相对容差。
updating, optional
如果是'immediate'
,则最佳解向量在单一生成过程中持续更新[4]。这可以加快收敛速度,因为试验向量可以利用最佳解的持续改进。如果是'deferred'
,则最佳解向量每代更新一次。只有'deferred'
与并行化或向量化兼容,且workers和vectorized关键字可以覆盖此选项。
自 1.2.0 版新增。
工作者int 或类似映射的可调用对象,可选
如果workers是整数,则将种群细分为workers部分,并并行评估(使用multiprocessing.Pool
)。提供-1 以使用所有可用 CPU 核心。或者提供类似映射的可调用对象,例如multiprocessing.Pool.map以并行评估种群。此评估以workers(func, iterable)
进行。如果workers != 1
,此选项将覆盖updating关键字为updating='deferred'
。如果workers != 1
,此选项将覆盖vectorized关键字。要求func可 pickle 化。
自 1.2.0 版新增。
约束条件
解算器的约束条件,除了bounds kwd 应用的约束外。使用 Lampinen 的方法[5]。
自 1.4.0 版新增。
x0None 或类似数组,可选
提供最小化的初始猜测。一旦种群被初始化,此向量替换第一个(最佳)成员。即使init给出了初始种群,也会进行此替换。x0.shape == (N,)
。
自 1.7.0 版新增。
整数性1-D 数组,可选
对于每个决策变量,一个布尔值,指示决策变量是否约束为整数值。该数组广播到(N,)
。如果有任何决策变量被约束为整数,则在精炼过程中它们不会改变。只使用介于下限和上限之间的整数值。如果在边界之间没有整数值,则会引发ValueError。
自 1.9.0 版新增。
向量化bool,可选
如果vectorized
为 True,则将使用形状为(N, S)
的x数组发送给func,并期望返回一个形状为(S,)
的数组,其中S是要计算的解向量数量。如果应用了约束条件,则用于构建Constraint对象的每个函数应接受形状为(N, S)
的x数组,并返回形状为(M, S)
的数组,其中M是约束组件的数量。这个选项是workers提供的并行化的替代方案,并且可以通过减少从多个函数调用中的解释器开销来提高优化速度。如果workers != 1
,则将忽略此关键字。此选项将覆盖updating关键字为updating='deferred'
。有关何时使用'vectorized'
和何时使用'workers'
的进一步讨论,请参见备注部分。
新版本 1.9.0 中的更新。
返回:
resOptimizeResult
优化结果表示为OptimizeResult
对象。重要属性包括:x
解向量数组,success
布尔标志,指示优化器是否成功退出,message
描述终止原因,population
种群中的解向量以及population_energies
每个population
条目的目标函数值。有关其他属性的描述,请参见OptimizeResult
。如果使用了polish,并且通过打磨获得了更低的最小值,则 OptimizeResult 还包含jac
属性。如果最终解决方案不满足应用的约束条件,则success
将为False。
备注
差分进化是一种基于随机种群的方法,适用于全局优化问题。在每次种群遍历过程中,算法通过与其他候选解混合来生成一个试验候选解,对每个候选解进行突变。有几种策略[3]用于创建试验候选解,适用于不同的问题。"best1bin"策略对许多系统而言是一个良好的起点。在这种策略中,随机选择种群中的两个成员。他们的差异用于突变迄今为止最佳成员("best"在"best1bin"中的意思),(x_0):
[b' = x_0 + mutation * (x_{r_0} - x_{r_1})]
然后构造一个试验向量。从随机选择的第 i 个参数开始,依次用来自b'
或原始候选者的参数填充(取模)。使用b'
或原始候选者的选择是通过二项分布('best1bin'中的'bin')进行的 - 生成一个介于[0, 1)的随机数。如果此数小于recombination常数,则参数从b'
加载,否则从原始候选者加载。最后一个参数始终从b'
加载。构建完试验候选者后,评估其适应度。如果试验结果比原始候选者好,则替换它。如果它还优于整体最佳候选者,则也替换该候选者。
其他可用的策略见 Qiang 和 Mitchell(2014)[3]。
[ \begin{align}\begin{aligned}rand1* : b' = x_{r_0} + mutation(x_{r_1} - x_{r_2})\rand2 : b' = x_{r_0} + mutation(x_{r_1} + x_{r_2} - x_{r_3} - x_{r_4})\best1 : b' = x_0 + mutation(x_{r_0} - x_{r_1})\best2 : b' = x_0 + mutation(x_{r_0} + x_{r_1} - x_{r_2} - x_{r_3})\currenttobest1 : b' = x_i + mutation(x_0 - x_i + x_{r_0} - x_{r_1})\randtobest1 : b' = x_{r_0} + mutation*(x_0 - x_{r_0} + x_{r_1} - x_{r_2})\end{aligned}\end{align} ]
其中整数(r_0, r_1, r_2, r_3, r_4)从区间[0, NP)随机选择,其中NP是总体大小,原始候选者索引为i。用户可以通过向strategy
提供可调用对象来完全自定义试验候选者的生成方式。
要提高找到全局最小值的机会,可以使用更高的popsize值,更高的mutation值和(抖动),但更低的recombination值。这样做会扩大搜索半径,但会减慢收敛速度。
默认情况下,最佳解向量在单次迭代中持续更新(updating='immediate'
)。这是对原始差分进化算法的修改[4],可以使试验向量立即从改进的解决方案中获益。要使用斯托恩和普莱斯的原始行为,即每次迭代更新一次最佳解,设置updating='deferred'
。'deferred'
方法既兼容并行化和矢量化('workers'
和'vectorized'
关键字)。这些方法可能通过更有效地利用计算资源来提高最小化速度。'workers'
将计算分布到多个处理器上。默认情况下,Python 的multiprocessing
模块用于此,但也可以使用其他方法,例如在集群上使用消息传递接口(MPI)[6] [7]。这些方法(创建新进程等)可能会带来一定的开销,这意味着计算速度并不一定随使用处理器数量的增加而线性扩展。并行化最适合计算密集型目标函数。如果目标函数较为简单,则'vectorized'
可能会有所帮助,因为它每次迭代仅调用一次目标函数,而不是为所有种群成员多次调用;解释器开销减少。
0.15.0 版新增内容。
参考文献
[1]
差分进化,维基百科,en.wikipedia.org/wiki/Differential_evolution
[2]
斯托恩(Storn, R)和普莱斯(Price, K),差分进化 - 一种简单高效的连续空间全局优化启发式算法,全球优化期刊,1997 年,11,341 - 359。
[3] (1,2)
强(Qiang, J.),米切尔(Mitchell, C.),统一差分进化算法用于全局优化,2014 年,www.osti.gov/servlets/purl/1163659
[4] (1,2)
沃明顿(Wormington, M.),帕纳乔内(Panaccione, C.),马特尼(Matney, K. M.),鲍文(Bowen, D. K.),利用遗传算法从 X 射线散射数据中表征结构,伦敦皇家学会 A 类学报,1999 年,357,2827-2848
[5]
兰皮宁(Lampinen, J.),差分进化算法的约束处理方法。2002 年进化计算大会论文集。CEC’02(Cat. No. 02TH8600)。第 2 卷。IEEE,2002 年。
[6]
mpi4py.readthedocs.io/en/stable/
[7]
schwimmbad.readthedocs.io/en/latest/
例子
让我们考虑最小化 Rosenbrock 函数的问题。这个函数在 rosen
中实现于 scipy.optimize
。
>>> import numpy as np
>>> from scipy.optimize import rosen, differential_evolution
>>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
>>> result = differential_evolution(rosen, bounds)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
现在重复,但使用并行化。
>>> result = differential_evolution(rosen, bounds, updating='deferred',
... workers=2)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)
让我们进行有约束的最小化。
>>> from scipy.optimize import LinearConstraint, Bounds
我们增加约束条件,即 x[0]
和 x[1]
的和必须小于或等于 1.9。这是一个线性约束,可以写成 A @ x <= 1.9
,其中 A = array([[1, 1]])
。这可以编码为 LinearConstraint
实例:
>>> lc = LinearConstraint([[1, 1]], -np.inf, 1.9)
使用 Bounds
对象指定限制。
>>> bounds = Bounds([0., 0.], [2., 2.])
>>> result = differential_evolution(rosen, bounds, constraints=lc,
... seed=1)
>>> result.x, result.fun
(array([0.96632622, 0.93367155]), 0.0011352416852625719)
寻找 Ackley 函数的最小值(en.wikipedia.org/wiki/Test_functions_for_optimization
)。
>>> def ackley(x):
... arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
... arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
... return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
>>> bounds = [(-5, 5), (-5, 5)]
>>> result = differential_evolution(ackley, bounds, seed=1)
>>> result.x, result.fun
(array([0., 0.]), 4.440892098500626e-16)
Ackley 函数以矢量化方式编写,因此可以使用 'vectorized'
关键字。请注意减少的函数评估次数。
>>> result = differential_evolution(
... ackley, bounds, vectorized=True, updating='deferred', seed=1
... )
>>> result.x, result.fun
(array([0., 0.]), 4.440892098500626e-16)
以下自定义策略函数模仿 'best1bin':
>>> def custom_strategy_fn(candidate, population, rng=None):
... parameter_count = population.shape(-1)
... mutation, recombination = 0.7, 0.9
... trial = np.copy(population[candidate])
... fill_point = rng.choice(parameter_count)
...
... pool = np.arange(len(population))
... rng.shuffle(pool)
...
... # two unique random numbers that aren't the same, and
... # aren't equal to candidate.
... idxs = []
... while len(idxs) < 2 and len(pool) > 0:
... idx = pool[0]
... pool = pool[1:]
... if idx != candidate:
... idxs.append(idx)
...
... r0, r1 = idxs[:2]
...
... bprime = (population[0] + mutation *
... (population[r0] - population[r1]))
...
... crossovers = rng.uniform(size=parameter_count)
... crossovers = crossovers < recombination
... crossovers[fill_point] = True
... trial = np.where(crossovers, bprime, trial)
... return trial
scipy.optimize.shgo
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.shgo.html#scipy.optimize.shgo
scipy.optimize.shgo(func, bounds, args=(), constraints=None, n=100, iters=1, callback=None, minimizer_kwargs=None, options=None, sampling_method='simplicial', *, workers=1)
使用 SHG 优化找到函数的全局最小值。
SHGO 代表“单纯同调全局优化”。
参数:
func可调用
要最小化的目标函数。必须以f(x, *args)
形式,其中x
是一个 1-D 数组的参数,而args
是完全指定函数所需的额外固定参数的元组。
bounds序列或Bounds
变量的边界。有两种指定边界的方法:
-
Bounds
类的实例。 -
x的每个元素的
(min, max)
对的序列。
args元组,可选
任何完全指定目标函数所需的额外固定参数。
constraints或约束列表,可选
约束定义。仅适用于 COBYLA、SLSQP 和 trust-constr。参见教程[5]以了解有关指定约束的详细信息。
注意
目前仅 COBYLA、SLSQP 和 trust-constr 局部最小化方法支持约束参数。如果在本地优化问题中使用的constraints
序列未在minimizer_kwargs
中定义,并且使用了约束方法,则将使用全局的constraints
(在minimizer_kwargs
中定义了constraints
序列意味着不会添加constraints
,因此如果需要添加等式约束等等,则需要将约束函数添加到minimizer_kwargs
中的不等式函数中)。COBYLA 仅支持不等式约束。
从版本 1.11.0 起:constraints
接受NonlinearConstraint
,LinearConstraint
。
n整数,可选
在构建单纯复合物时使用的采样点数。对于默认的simplicial
采样方法,生成 2**dim + 1 个采样点,而不是默认的n=100。对于所有其他指定的值,生成n个采样点。对于sobol
、halton
和其他任意的sampling_methods,生成n=100或另一个指定的采样点数。
iters整数,可选
用于构建单纯复合物的迭代次数。默认为 1。
callback可调用,可选
在每次迭代后调用,形式为callback(xk)
,其中xk
是当前参数向量。
minimizer_kwargs字典,可选
要传递给最小化器scipy.optimize.minimize
的额外关键字参数。一些重要的选项可能是:
methodstr
最小化方法。如果未指定,则根据问题是否有约束或边界选择为 BFGS、L-BFGS-B、SLSQP 之一。
argstuple
传递给目标函数 (
func
) 及其导数(Jacobian、Hessian)的额外参数。选项字典,可选
注意,默认情况下容差被指定为
{ftol: 1e-12}
选项字典,可选
解算器选项字典。许多指定给全局例程的选项也传递给 scipy.optimize.minimize 例程。传递给局部例程的选项标有“(L)”。
停止标准,如果满足任何指定的准则则算法将终止。但是,默认算法不需要指定任何准则:
-
maxfevint (L)
可行域内的最大函数评估次数。(注意,只有支持此选项的方法才会精确地在指定值处终止例程。否则,准则只会在全局迭代期间终止)
-
f_min
如果已知,指定最小目标函数值。
-
f_tolfloat
值 f 的停止准则的精度目标。请注意,如果全局例程中的采样点在此容差内,则全局例程也将终止。
-
maxiterint
执行的最大迭代次数。
-
maxevint
最大采样评估次数(包括在不可行点中的搜索)。
-
maxtimefloat
允许的最大处理运行时
-
minhgrdint
最小同调群秩差分。在每次迭代期间(大约)计算目标函数的同调群。该群的秩与目标函数中局部凸子域的数量具有一一对应关系(在足够的采样点后,这些子域包含唯一的全局最小值)。如果指定迭代的
maxhgrd
中 hgr 的差异为 0,则算法将终止。
目标函数知识:
-
对称性列表或布尔值
指定目标函数是否包含对称变量。在完全对称的情况下,搜索空间(因此性能)最多减少 O(n!) 倍。如果指定为 True,则所有变量将被设置为相对于第一个变量对称。默认设置为 False。
例如,f(x) = (x_1 + x_2 + x_3) + (x_4)2 + (x_5)2 + (x_6)**2
在此方程中,x_2 和 x_3 对 x_1 对称,而 x_5 和 x_6 对 x_4 对称,可以指定给解算器如下:
对称性 = [0, # 变量 1
0, # 对变量 1 对称 0, # 对变量 1 3, # 变量 4 3, # 对变量 4 3, # 对变量 4 ]
-
jacbool 或可调用,可选
目标函数的 Jacobian(梯度)。仅适用于 CG、BFGS、Newton-CG、L-BFGS-B、TNC、SLSQP、dogleg 和 trust-ncg。如果
jac
为布尔值且为 True,则假定fun
返回目标函数的梯度。如果为 False,则梯度将以数值方式估计。jac
也可以是一个返回目标函数梯度的可调用对象。在这种情况下,它必须接受与fun
相同的参数。(将自动传递给scipy.optimize.minimize
) -
hess,hesspcallable,可选
Hessian(二阶导数矩阵)的目标函数或者目标函数的 Hessian 乘以任意向量 p。仅适用于 Newton-CG、dogleg 和 trust-ncg。
hessp
或hess
中只需提供一个。如果提供了hess
,则将忽略hessp
。如果hess
和hessp
都未提供,则将在jac
上使用有限差分近似 Hessian 乘积。hessp
必须计算 Hessian 乘以任意向量。(将自动传递给scipy.optimize.minimize
)
算法设置:
-
minimize_every_iterbool
如果为 True,则有前景的全局采样点将传递给每次迭代的本地最小化程序。如果为 False,则仅运行最终的最小化池。默认为 True。
-
local_iterint
每次迭代仅评估少数最佳最小化候选点。如果为 False,则所有潜在点都将传递给本地最小化程序。
-
infty_constraintsbool
如果为 True,则生成的任何采样点超出可行域将被保存,并给予目标函数值
inf
。如果为 False,则这些点将被丢弃。使用此功能可以在找到全局最小值之前提高函数评估的性能。指定为 False 将以稍微降低性能的代价节省内存。默认为 True。
反馈:
-
dispbool(L)
设置为 True 以打印收敛消息。
sampling_methodstr 或函数,可选
当前内置的采样方法选项为halton
、sobol
和simplicial
。默认的simplicial
提供了有限时间内收敛到全局最小值的理论保证。halton
和sobol
方法在采样点生成方面更快,但失去了保证收敛性。对于大多数“较简单”的问题,这更为适用。用户定义的采样函数必须每次调用接受n
个维度为dim
的采样点,并输出形状为n x dim的采样点数组。
workersint 或类似映射的可调用对象,可选
并行地对样本进行本地串行最小化。提供-1 以使用所有可用的 CPU 核心,或提供一个整数以使用这么多进程(使用multiprocessing.Pool
)。
或者提供一个类似映射的可调用对象,例如multiprocessing.Pool.map以进行并行评估。此评估以workers(func, iterable)
形式进行。要求func可被 pickle。
从版本 1.11.0 开始新增。
返回:
resOptimizeResult
优化结果表示为OptimizeResult
对象。重要属性包括:x
是对应全局最小值的解数组,fun
是全局解处的函数输出,xl
是局部最小解的排序列表,funl
是相应局部解的函数输出,success
是一个布尔标志,指示优化器是否成功退出,message
描述终止原因,nfev
是包括采样调用在内的总目标函数评估次数,nlfev
是所有局部搜索优化导致的总目标函数评估次数,nit
是全局例程执行的迭代次数。
注释
使用单纯同调全局优化进行全局优化 [1]。适用于解决通用 NLP 和黑盒优化问题以达到全局最优(低维问题)。
一般来说,优化问题的形式为:
minimize f(x) subject to
g_i(x) >= 0, i = 1,...,m
h_j(x) = 0, j = 1,...,p
这里 x 是一个或多个变量的向量。f(x)
是目标函数R^n -> R
,g_i(x)
是不等式约束,h_j(x)
是等式约束。
可选地,还可以使用bounds参数指定 x 中每个元素的下限和上限。
虽然 SHGO 的大部分理论优势仅对f(x)
为 Lipschitz 光滑函数时得到证明,但是当f(x)
是非连续、非凸且非光滑时,如果使用默认的采样方法,该算法也被证明能够收敛到全局最优解 [1]。
可以使用minimizer_kwargs
参数指定本地搜索方法,该参数传递给scipy.optimize.minimize
。默认情况下使用SLSQP
方法。一般建议如果问题定义了不等式约束,则使用SLSQP
或COBYLA
本地最小化方法,因为其他方法不使用约束。
使用scipy.stats.qmc
生成halton
和sobol
方法点。还可以使用任何其他 QMC 方法。
参考文献
[1] (1,2)
Endres, SC, Sandrock, C, Focke, WW (2018) “一种用于 Lipschitz 优化的单纯同调算法”,全球优化期刊。
[2]
Joe, SW 和 Kuo, FY(2008)“用更好的二维投影构建 Sobol 序列”,SIAM J. Sci. Comput. 30, 2635-2654。
[3] (1,2)
Hock, W 和 Schittkowski, K(1981)“非线性规划代码的测试示例”,经济与数学系统讲义,187. Springer-Verlag,纽约。www.ai7.uni-bayreuth.de/test_problem_coll.pdf
[4]
Wales, DJ(2015)“观点:从势能景观中获取反应坐标和动态的洞察”,化学物理学杂志,142(13), 2015。
[5]
示例
首先考虑最小化 Rosenbrock 函数的问题,rosen
:
>>> from scipy.optimize import rosen, shgo
>>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
>>> result = shgo(rosen, bounds)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 2.920392374190081e-18)
注意,边界确定目标函数的维度,因此是必需的输入,但是您可以使用None
或类似np.inf
的对象指定空边界,这些将被转换为大的浮点数。
>>> bounds = [(None, None), ]*4
>>> result = shgo(rosen, bounds)
>>> result.x
array([0.99999851, 0.99999704, 0.99999411, 0.9999882 ])
接下来,我们考虑 Eggholder 函数,这是一个具有多个局部极小值和一个全局极小值的问题。我们将演示shgo
的参数使用和能力。
>>> import numpy as np
>>> def eggholder(x):
... return (-(x[1] + 47.0)
... * np.sin(np.sqrt(abs(x[0]/2.0 + (x[1] + 47.0))))
... - x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47.0))))
... )
...
>>> bounds = [(-512, 512), (-512, 512)]
shgo
具有内置的低差异采样序列。首先,我们将输入Sobol'序列的 64 个初始采样点:
>>> result = shgo(eggholder, bounds, n=64, sampling_method='sobol')
>>> result.x, result.fun
(array([512\. , 404.23180824]), -959.6406627208397)
shgo
还返回了找到的任何其他局部极小值,可以使用以下方式调用:
>>> result.xl
array([[ 512\. , 404.23180824],
[ 283.0759062 , -487.12565635],
[-294.66820039, -462.01964031],
[-105.87688911, 423.15323845],
[-242.97926 , 274.38030925],
[-506.25823477, 6.3131022 ],
[-408.71980731, -156.10116949],
[ 150.23207937, 301.31376595],
[ 91.00920901, -391.283763 ],
[ 202.89662724, -269.38043241],
[ 361.66623976, -106.96493868],
[-219.40612786, -244.06020508]])
>>> result.funl
array([-959.64066272, -718.16745962, -704.80659592, -565.99778097,
-559.78685655, -557.36868733, -507.87385942, -493.9605115 ,
-426.48799655, -421.15571437, -419.31194957, -410.98477763])
这些结果在应用中非常有用,特别是在需要许多全局极小值和其他全局极小值值的情况下,或者在局部极小值可以为系统提供洞察力的情况下(例如物理化学中的形态学[4])。
如果我们想要找到更多的局部极小值,我们可以增加采样点或迭代次数的数量。我们将增加采样点数到 64,并将迭代次数从默认值 1 增加到 3。使用simplicial
,这将为我们提供 64 x 3 = 192 个初始采样点。
>>> result_2 = shgo(eggholder,
... bounds, n=64, iters=3, sampling_method='sobol')
>>> len(result.xl), len(result_2.xl)
(12, 23)
注意,例如n=192, iters=1
和n=64, iters=3
之间的差异。在第一种情况下,仅一次处理最小化池中的有前途点。在后一种情况下,它每 64 个采样点处理一次,总共 3 次。
要演示解决具有非线性约束的问题,请考虑 Hock 和 Schittkowski 问题 73(牛饲料)的以下示例[3]:
minimize: f = 24.55 * x_1 + 26.75 * x_2 + 39 * x_3 + 40.50 * x_4
subject to: 2.3 * x_1 + 5.6 * x_2 + 11.1 * x_3 + 1.3 * x_4 - 5 >= 0,
12 * x_1 + 11.9 * x_2 + 41.8 * x_3 + 52.1 * x_4 - 21
-1.645 * sqrt(0.28 * x_1**2 + 0.19 * x_2**2 +
20.5 * x_3**2 + 0.62 * x_4**2) >= 0,
x_1 + x_2 + x_3 + x_4 - 1 == 0,
1 >= x_i >= 0 for all i
在[3]中给出的近似答案是:
f([0.6355216, -0.12e-11, 0.3127019, 0.05177655]) = 29.894378
>>> def f(x): # (cattle-feed)
... return 24.55*x[0] + 26.75*x[1] + 39*x[2] + 40.50*x[3]
...
>>> def g1(x):
... return 2.3*x[0] + 5.6*x[1] + 11.1*x[2] + 1.3*x[3] - 5 # >=0
...
>>> def g2(x):
... return (12*x[0] + 11.9*x[1] +41.8*x[2] + 52.1*x[3] - 21
... - 1.645 * np.sqrt(0.28*x[0]**2 + 0.19*x[1]**2
... + 20.5*x[2]**2 + 0.62*x[3]**2)
... ) # >=0
...
>>> def h1(x):
... return x[0] + x[1] + x[2] + x[3] - 1 # == 0
...
>>> cons = ({'type': 'ineq', 'fun': g1},
... {'type': 'ineq', 'fun': g2},
... {'type': 'eq', 'fun': h1})
>>> bounds = [(0, 1.0),]*4
>>> res = shgo(f, bounds, n=150, constraints=cons)
>>> res
message: Optimization terminated successfully.
success: True
fun: 29.894378159142136
funl: [ 2.989e+01]
x: [ 6.355e-01 1.137e-13 3.127e-01 5.178e-02] # may vary
xl: [[ 6.355e-01 1.137e-13 3.127e-01 5.178e-02]] # may vary
nit: 1
nfev: 142 # may vary
nlfev: 35 # may vary
nljev: 5
nlhev: 0
>>> g1(res.x), g2(res.x), h1(res.x)
(-5.062616992290714e-14, -2.9594104944408173e-12, 0.0)
scipy.optimize.dual_annealing
scipy.optimize.dual_annealing(func, bounds, args=(), maxiter=1000, minimizer_kwargs=None, initial_temp=5230.0, restart_temp_ratio=2e-05, visit=2.62, accept=-5.0, maxfun=10000000.0, seed=None, no_local_search=False, callback=None, x0=None)
使用双重退火法找到函数的全局最小值。
参数:
func可调用对象
要最小化的目标函数。必须以f(x, *args)
的形式给出,其中x
是一维数组形式的参数,args
是一个包含完全指定函数所需的任何额外固定参数的元组。
边界序列或Bounds
类
变量的边界。有两种指定边界的方式:
-
Bounds
类的实例。 -
对于x中的每个元素,都有
(min, max)
对的序列。
args元组, 可选
任何完全指定目标函数所需的额外固定参数。
最大迭代次数int, 可选
全局搜索迭代的最大次数。默认值为 1000。
minimizer_kwargs字典, 可选
传递给局部最小化器(minimize
)的额外关键字参数。一些重要的选项可能包括:method
用于指定最小化方法和args
用于目标函数的额外参数。
初始温度float, 可选
初始温度,使用较高的值可以促进更广泛的能量景观搜索,允许dual_annealing
逃离被困在其中的局部极小值。默认值为 5230。范围为(0.01, 5.e4]。
重启温度比率float, 可选
在退火过程中,温度逐渐降低,当达到initial_temp * restart_temp_ratio
时,会触发重新退火过程。比率的默认值为 2e-5。范围为(0, 1)。
访问率float, 可选
访问分布参数。默认值为 2.62。较高的值使访问分布尾部更重,这使得算法跳到更远的区域。值的范围为(1, 3]。
接受率float, 可选
接受分布参数。用于控制接受概率。接受参数越低,接受概率越小。默认值为-5.0,范围为(-1e4, -5]。
最大函数调用次数int, 可选
目标函数调用的软限制。如果算法在局部搜索中间,超出这个数值后,算法将在局部搜索完成后停止。默认值为 1e7。
种子{None, int, numpy.random.Generator
, numpy.random.RandomState
}, 可选
如果 seed 为 None(或 np.random),则使用 numpy.random.RandomState
单例。如果 seed 是一个整数,则使用一个新的 RandomState
实例,并用 seed 初始化。如果 seed 已经是 Generator
或 RandomState
实例,则使用该实例。指定 seed 可重复进行最小化。使用此种子生成的随机数仅影响访问分布函数和新坐标生成。
no_local_searchbool, optional
如果将 no_local_search 设置为 True,则将执行传统的广义模拟退火,不会应用局部搜索策略。
callbackcallable, optional
带有签名 callback(x, f, context)
的回调函数,将对找到的所有最小值进行调用。x
和 f
是最新最小值的坐标和函数值,context
的值在 [0, 1, 2] 范围内,具有以下含义:
- 0:在退火过程中检测到最小值。
- 1:在局部搜索过程中发生检测。
- 2:在双退火过程中完成检测。
如果回调实现返回 True,则算法将停止。
x0ndarray, shape(n,), optional
单个 N-D 起始点的坐标。
返回:
resOptimizeResult
优化结果表示为 OptimizeResult
对象。重要属性包括:x
解数组,fun
在解处的函数值,以及 message
描述终止原因。查看 OptimizeResult
以了解其他属性的描述。
注:
这个函数实现了双退火优化。这种随机方法源于[3],结合了 CSA(经典模拟退火)和 FSA(快速模拟退火)[1] [2],并采用一种策略,在接受的位置上应用局部搜索[4]。这种算法的另一种实现在[5]中有描述,并在[6]中进行了基准测试。这种方法引入了一种高级方法来优化广义退火过程中找到的解。该算法使用扭曲的柯西-洛伦兹访问分布,其形状由参数 (q_{v}) 控制。
[g_{q_{v}}(\Delta x(t)) \propto \frac{ \ \left[T_{q_{v}}(t) \right]^{-\frac{D}{3-q_{v}}}}{ \ \left[{1+(q_{v}-1)\frac{(\Delta x(t))^{2}} { \ \left[T_{q_{v}}(t)\right]{\frac{2}{3-q_{v}}}}}\right]-1}+\frac{D-1}{2}}}]
其中 (t) 是人造时间。这种访问分布用于生成变量 (x(t)) 的试验跳跃距离 (\Delta x(t)),在人造温度 (T_{q_{v}}(t)) 下。
从起始点开始,调用访问分布函数后,接受概率计算如下:
[p_{q_{a}} = \min{{1,\left[1-(1-q_{a}) \beta \Delta E \right]^{ \ \frac{1}{1-q_{a}}}}}]
其中 (q_{a}) 是接受参数。对于 (q_{a}<1),在 (1-(1-q_{a}) \beta \Delta E < 0) 的情况下,将分配零接受概率。
[[1-(1-q_{a}) \beta \Delta E] < 0]
人工温度 (T_{q_{v}}(t)) 根据以下方式递减:
[T_{q_{v}}(t) = T_{q_{v}}(1) \frac{2^{q_{v}-1}-1}{\left( \ 1 + t\right)^{q_{v}-1}-1}]
其中 (q_{v}) 是访问参数。
版本 1.2.0 的新功能。
参考文献
[1]
Tsallis C. 可能是 Boltzmann-Gibbs 统计的一般化。《统计物理学杂志》,52, 479-487 (1998)。
[2]
Tsallis C, Stariolo DA. 广义模拟退火。《物理学 A》,233, 395-406 (1996)。
[3]
Xiang Y, Sun DY, Fan W, Gong XG. 广义模拟退火算法及其在汤姆森模型中的应用。《物理学快报 A》,233, 216-220 (1997)。DOI:10.18637/jss.v060.i06
[4]
Xiang Y, Gong XG. 广义模拟退火算法的效率。《物理评论 E》,62, 4473 (2000)。
[5]
Xiang Y, Gubian S, Suomela B, Hoeng J. 用于高效全局优化的广义模拟退火:R 包 GenSA。《R 语言杂志》,Volume 5/1 (2013)。
[6]
Mullen, K. R 中的连续全局优化。《统计软件杂志》,60(6), 1 - 45, (2014)。DOI:10.18637/jss.v060.i06
示例
下面的例子是一个 10 维问题,有许多局部极小值。涉及的函数称为 Rastrigin (en.wikipedia.org/wiki/Rastrigin_function
)
>>> import numpy as np
>>> from scipy.optimize import dual_annealing
>>> func = lambda x: np.sum(x*x - 10*np.cos(2*np.pi*x)) + 10*np.size(x)
>>> lw = [-5.12] * 10
>>> up = [5.12] * 10
>>> ret = dual_annealing(func, bounds=list(zip(lw, up)))
>>> ret.x
array([-4.26437714e-09, -3.91699361e-09, -1.86149218e-09, -3.97165720e-09,
-6.29151648e-09, -6.53145322e-09, -3.93616815e-09, -6.55623025e-09,
-6.05775280e-09, -5.00668935e-09]) # random
>>> ret.fun
0.000000
scipy.optimize.direct
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.direct.html#scipy.optimize.direct
scipy.optimize.direct(func, bounds, *, args=(), eps=0.0001, maxfun=None, maxiter=1000, locally_biased=True, f_min=-inf, f_min_rtol=0.0001, vol_tol=1e-16, len_tol=1e-06, callback=None)
使用 DIRECT 算法找到函数的全局最小值。
参数:
funccallable
要最小化的目标函数。func(x, *args) -> float
,其中x
是形状为(n,)的一维数组,args
是完全指定函数所需的固定参数的元组。
bounds序列或Bounds
变量的界限。有两种指定界限的方式:
-
Bounds
类的实例。 -
对于
x
中的每个元素,(min, max)
对。
argstuple,可选
完全指定目标函数所需的任何额外固定参数。
epsfloat,可选
当前最佳超矩形与下一个可能的最优超矩形之间的目标函数值的最小必需差异。因此,eps用作局部和全局搜索之间的折衷:它越小,搜索就越局部。默认为 1e-4。
maxfunint 或 None,可选
目标函数评估的大致上限。如果为None,则自动设置为1000 * N
,其中N
表示维度的数量。如有必要,将对 DIRECT 的 RAM 使用进行限制,以保持约为 1GiB。这仅适用于维度非常高和max_fun过大的问题。默认为None。
maxiterint,可选
最大迭代次数。默认为 1000。
locally_biasedbool,可选
如果为True(默认值),则使用称为 DIRECT_L 的算法的局部偏置变体。如果为False,则使用原始的无偏 DIRECT 算法。对于具有许多局部最小值的困难问题,建议使用False。
f_minfloat,可选
全局最优解的函数值。仅在已知全局最优解时设置此值。默认为-np.inf
,因此此终止准则被禁用。
f_min_rtolfloat,可选
当当前最佳最小值f和提供的全局最小值f_min之间的相对误差小于f_min_rtol时,终止优化。如果f_min也设置了,则使用此参数。必须介于 0 和 1 之间。默认为 1e-4。
vol_tolfloat,可选
当包含最低函数值的超矩形的体积小于完整搜索空间的vol_tol时,终止优化。必须介于 0 和 1 之间。默认为 1e-16。
len_tolfloat,可选
如果 locally_biased=True,则当包含最低函数值的超矩形的归一化最大边长的一半小于 len_tol 时终止优化。如果 locally_biased=False,则当包含最低函数值的超矩形的归一化对角线的一半小于 len_tol 时终止优化。必须介于 0 和 1 之间。默认值为 1e-6。
callback 可调用对象,可选
具有签名 callback(xk)
的回调函数,其中 xk
表示迄今为止找到的最佳函数值。
返回:
res OptimizeResult
表示优化结果的 OptimizeResult
对象。重要属性包括:x
解数组,success
布尔标志,指示优化器是否成功退出,以及 message
描述终止原因。详见 OptimizeResult
获取其他属性的描述。
笔记
DIviding RECTangles (DIRECT) 是一种确定性全局优化算法,能够通过在搜索空间中采样潜在解来最小化黑盒函数,其中变量受下限和上限约束 [1]。该算法首先将搜索空间标准化为 n 维单位超立方体。它在这个超立方体的中心点和每个坐标方向上的 2n 个点处采样函数。使用这些函数值,DIRECT 将域划分为超矩形,每个超矩形的中心点恰好是一个采样点。在每次迭代中,DIRECT 使用默认为 1e-4 的eps参数选择一些现有超矩形进行进一步划分。这个划分过程持续进行,直到达到最大迭代次数或允许的最大函数评估次数,或者包含到目前为止找到的最小值的超矩形足够小。如果指定了f_min,优化将在相对容差内达到这个函数值时停止。默认情况下使用 DIRECT 的局部偏向变体(最初称为 DIRECT_L) [2]。它使搜索更加局部偏向,并且对只有少数局部最小值的情况更有效。
关于终止标准的说明:vol_tol 是指包含到目前为止找到的最低函数值的超矩形的体积。这个体积随问题维数的增加呈指数下降。因此,为了避免在更高维问题上过早终止算法,应减小 vol_tol。但对于 len_tol 不适用此规则:它指的是最大边长的一半(对于 locally_biased=True
)或者超矩形对角线的一半(对于 locally_biased=False
)。
这段代码基于 Gablonsky 等人的 DIRECT 2.0.4 Fortran 代码,可以在ctk.math.ncsu.edu/SOFTWARE/DIRECTv204.tar.gz
找到。这个原始版本最初通过 f2c 转换,然后由 Steven G. Johnson 在 2007 年 8 月为 NLopt 项目进行了清理和重新组织。direct
函数封装了 C 实现。
新功能版本 1.9.0 中的更新。
参考文献
[1]
Jones, D.R., Perttunen, C.D. & Stuckman, B.E. 没有 Lipschitz 常数的 Lipschitz 优化。《优化理论与应用期刊》79, 157-181 (1993).
[2]
Gablonsky, J., Kelley, C. DIRECT 算法的一种本地偏向形式。《全局优化期刊》21, 27-37 (2001).
示例
以下示例是一个二维问题,有四个局部最小值:最小化 Styblinski-Tang 函数(en.wikipedia.org/wiki/Test_functions_for_optimization
)。
>>> from scipy.optimize import direct, Bounds
>>> def styblinski_tang(pos):
... x, y = pos
... return 0.5 * (x**4 - 16*x**2 + 5*x + y**4 - 16*y**2 + 5*y)
>>> bounds = Bounds([-4., -4.], [4., 4.])
>>> result = direct(styblinski_tang, bounds)
>>> result.x, result.fun, result.nfev
array([-2.90321597, -2.90321597]), -78.3323279095383, 2011
找到了正确的全局最小值,但使用了大量的函数评估(2011)。可以通过放宽终止容差vol_tol和len_tol来提前停止 DIRECT。
>>> result = direct(styblinski_tang, bounds, len_tol=1e-3)
>>> result.x, result.fun, result.nfev
array([-2.9044353, -2.9044353]), -78.33230330754142, 207
scipy.optimize.least_squares
scipy.optimize.least_squares(fun, x0, jac='2-point', bounds=(-inf, inf), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options={}, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs={})
解决带有变量边界的非线性最小二乘问题。
给定残差 f(x)(一个 m-D 实函数的 n 个实变量)和损失函数 rho(s)(一个标量函数),least_squares
找到成本函数 F(x)的局部最小值:
minimize F(x) = 0.5 * sum(rho(f_i(x)**2), i = 0, ..., m - 1)
subject to lb <= x <= ub
损失函数 rho(s)的目的是减少离群值对解决方案的影响。
参数:
fun可调用
计算残差向量的函数,具有签名fun(x, *args, **kwargs)
,即最小化将相对于其第一个参数进行。传递给此函数的参数x
是形状为(n,)的 ndarray(即使对于 n=1,也不是标量)。它必须分配并返回形状为(m,)的 1-D 数组或标量。如果参数x
是复数或函数fun
返回复数残差,则必须将其包装在实数参数的实函数中,如示例部分的最后所示。
x0array_like,形状为(n,)或 float
独立变量的初始猜测。如果是 float,它将被视为具有一个元素的 1-D 数组。当method为‘trf’时,初始猜测可能会被略微调整,以确保足够位于给定的bounds内。
jac,可选
计算雅可比矩阵的方法(一个 m×n 矩阵,其中元素(i, j)是 f[i]关于 x[j]的偏导数)。关键字选择了用于数值估计的有限差分方案。方案“3-point”更精确,但需要两倍的操作量比“2-point”(默认)多。方案“cs”使用复步长,虽然可能是最精确的,但仅当fun正确处理复数输入并且能在复平面上解析延续时适用。方法“lm”始终使用“2-point”方案。如果可调用,则用作jac(x, *args, **kwargs)
,应返回雅可比矩阵的良好近似值(或确切值),作为数组(np.atleast_2d 应用),稀疏矩阵(出于性能考虑,优先使用 csr_matrix)或scipy.sparse.linalg.LinearOperator
。
bounds形如 2 元组的 array_like 或Bounds
,可选
有两种指定边界的方法:
Bounds
类的实例- 独立变量的上下界。默认无界。每个数组必须与x0的大小相匹配或者是标量,在后一种情况下,所有变量的边界将相同。使用适当符号的
np.inf
以禁用所有或部分变量的边界。
method, optional
执行最小化的算法。
- ‘trf’:反射信赖区域算法,特别适用于具有边界的大型稀疏问题。通常是健壮的方法。
- ‘dogbox’:矩形信赖区域的狗腿算法,典型应用场景是具有边界的小问题。不推荐处理秩亏的雅可比矩阵问题。
- ‘lm’:MINPACK 中实现的 Levenberg-Marquardt 算法。不处理边界和稀疏雅可比矩阵。通常是小型无约束问题的最有效方法。
默认为‘trf’。更多信息请参见注释。
ftolfloat 或 None,可选
终止由成本函数变化的公差。默认为 1e-8。当dF < ftol * F
且在上一步中局部二次模型与真实模型之间有足够协议时,优化过程停止。
若为 None 且‘method’不为‘lm’,则此条件下的终止被禁用。若‘method’为‘lm’,则此公差必须高于机器 epsilon。
xtolfloat 或 None,可选
终止由独立变量变化的公差。默认为 1e-8。确切条件取决于所使用的method:
- 对于‘trf’和‘dogbox’:
norm(dx) < xtol * (xtol + norm(x))
。- 对于‘lm’:
Delta < xtol * norm(xs)
,其中Delta
是信赖区域半径,xs
是根据x_scale参数缩放后的x的值(见下文)。
若为 None 且‘method’不为‘lm’,则此条件下的终止被禁用。若‘method’为‘lm’,则此公差必须高于机器 epsilon。
gtolfloat 或 None,可选
终止由梯度范数的公差。默认为 1e-8。确切条件取决于所使用的method:
- 对于‘trf’:
norm(g_scaled, ord=np.inf) < gtol
,其中g_scaled
是考虑边界存在的梯度值[STIR]。- 对于‘dogbox’:
norm(g_free, ord=np.inf) < gtol
,其中g_free
是相对于在边界上不处于最佳状态的变量的梯度。- 对于‘lm’:雅可比矩阵的列与残差向量之间夹角的最大绝对值小于gtol,或者残差向量为零。
若为 None 且‘method’不为‘lm’,则此条件下的终止被禁用。若‘method’为‘lm’,则此公差必须高于机器 epsilon。
x_scalearray_like 或‘jac’,可选
每个变量的特征尺度。设置 x_scale 相当于在缩放变量 xs = x / x_scale
中重新表述问题。另一种观点是,沿第 j 维的信任区域大小与 x_scale[j]
成比例。通过设置 x_scale,可以实现改进的收敛性,使得沿着任何缩放变量的给定大小步骤对成本函数产生类似的影响。如果设置为 ‘jac’,则使用雅可比矩阵列的逆范数进行迭代更新尺度(如 [JJMore] 中描述)。
lossstr 或 callable,可选
确定损失函数。允许以下关键字值:
- ‘linear’(默认):
rho(z) = z
。给出一个标准的最小二乘问题。- ‘soft_l1’:
rho(z) = 2 * ((1 + z)**0.5 - 1)
。l1(绝对值)损失的平滑逼近。通常是鲁棒最小二乘的良好选择。- ‘huber’ :
rho(z) = z if z <= 1 else 2*z**0.5 - 1
。与 ‘soft_l1’ 类似工作。- ‘cauchy’ :
rho(z) = ln(1 + z)
. 严重削弱离群值的影响,但可能会导致优化过程中的困难。- ‘arctan’:
rho(z) = arctan(z)
。限制单个残差的最大损失,具有类似 ‘cauchy’ 的特性。
如果可调用,它必须接受一个一维 ndarray z=f**2
并返回一个形状为 (3, m) 的 array_like,其中第 0 行包含函数值,第 1 行包含一阶导数,第 2 行包含二阶导数。方法 ‘lm’ 仅支持 ‘linear’ 损失。
f_scalefloat,可选
变量之间的软边界值,默认为 1.0。损失函数的评估如下 rho_(f**2) = C**2 * rho(f**2 / C**2)
,其中 C
是 f_scale,而 rho
受 loss 参数控制。对于 loss='linear'
,此参数没有影响,但对于其他 loss 值,它至关重要。
max_nfevNone 或 int,可选
函数评估的最大次数在终止之前。如果为 None(默认),则自动选择该值:
- 对于 ‘trf’ 和 ‘dogbox’:100 * n。
- For ‘lm’ : 100 * n if jac is callable and 100 * n * (n + 1) otherwise (because ‘lm’ counts function calls in Jacobian estimation).
diff_stepNone 或 array_like,可选
确定雅可比矩阵有限差分近似的相对步长。实际步长计算为 x * diff_step
。如果为 None(默认),则 diff_step 被认为是用于有限差分方案的常规“最优”机器 epsilon 的幂 [NR]。
tr_solver,可选
解决信任区域子问题的方法,仅对 ‘trf’ 和 ‘dogbox’ 方法有效。
- ‘exact’ 适用于问题规模不是很大且雅可比矩阵稠密的情况。每次迭代的计算复杂度与雅可比矩阵的奇异值分解相当。
- ‘lsmr’适用于具有稀疏和大雅可比矩阵的问题。它使用迭代过程
scipy.sparse.linalg.lsmr
来找到线性最小二乘问题的解,并且只需矩阵-向量乘积评估。
如果为 None(默认值),则求解器基于第一次迭代返回的雅可比矩阵类型选择。
tr_options字典,可选
关键字选项传递给信任区域求解器。
tr_solver='exact'
:tr_options将被忽略。tr_solver='lsmr'
:scipy.sparse.linalg.lsmr
的选项。此外,method='trf'
支持‘regularize’选项(布尔值,默认为 True),该选项在正规方程中添加正则化项,如果雅可比矩阵是秩缺乏的[Byrd](eq. 3.4)时可以改善收敛性。
jac_sparsity,可选
定义雅可比矩阵的稀疏结构用于有限差分估计,其形状必须为(m, n)。如果雅可比矩阵每行只有少量非零元素,则提供稀疏结构将极大加速计算[Curtis]。零条目意味着雅可比矩阵中对应元素恒为零。如果提供,则强制使用‘lsmr’信任区域求解器。如果为 None(默认值),则将使用密集差分。对‘lm’方法无效。
详细,可选
算法详细程度:
- 0(默认值):静默工作。
- 1:显示终止报告。
- 2:在迭代过程中显示进展(不支持‘lm’方法)。
args, kwargs元组和字典,可选
传递给fun和jac的额外参数。默认情况下两者为空。调用签名为fun(x, *args, **kwargs)
,jac相同。
返回:
resultOptimizeResult
OptimizeResult
,定义了以下字段:
xndarray,形状为(n,)
找到解决方案。
costfloat
解处成本函数的值。
funndarray,形状为(m,)
解处的残差向量。
jacndarray、稀疏矩阵或线性算子,形状为(m, n)
在解处的修改雅可比矩阵,即 J^T J 是成本函数海森矩阵的高斯-牛顿近似。类型与算法使用的类型相同。
gradndarray,形状为(m,)
解处成本函数的梯度。
优越性浮点数
一阶最优性度量。在无约束问题中,它总是梯度的一致范数。在有约束问题中,它是在迭代过程中与gtol比较的量。
active_maskndarray 整数,形状为(n,)
每个分量显示相应约束是否活跃(即变量是否位于边界处):
- 0:约束未生效。
- -1:下界生效。
- 1 : 上限激活。
对于 ‘trf’ 方法来说可能有些随意,因为它生成严格可行迭代序列,并且 active_mask 在容差阈值内确定。
nfevint
函数评估次数。‘trf’ 和 ‘dogbox’ 方法不计算数值雅可比逼近的函数调用次数,而‘lm’方法则计算。
njevint 或 None
雅可比矩阵评估次数。如果在 ‘lm’ 方法中使用数值雅可比逼近,则设置为 None。
statusint
算法终止的原因:
- -1 : 从 MINPACK 返回的不合适的输入参数状态。
- 0 : 超过了最大函数评估次数。
- 1 : gtol 终止条件满足。
- 2 : ftol 终止条件满足。
- 3 : xtol 终止条件满足。
- 4 : ftol 和 xtol 终止条件均满足。
messagestr
终止原因的口头描述。
successbool
如果满足其中一个收敛条件,则返回 True(status > 0)。
另见
leastsq
MINPACK 实现的 Levenberg-Marquadt 算法的遗留包装器。
curve_fit
用于曲线拟合问题的最小二乘最小化。
注释
‘lm’ 方法(Levenberg-Marquardt)调用 MINPACK 中实现的最小二乘算法的包装器(lmder, lmdif)。它以信任区域类型算法的形式运行 Levenberg-Marquardt 算法。该实现基于文献 [JJMore],非常稳健和高效,采用了许多巧妙的技巧。对于无约束问题,它应该是你的首选。注意它不支持边界条件,且在 m < n 时不起作用。
方法‘trf’(Trust Region Reflective)的动机源于解决方程组的过程,这组成了形式化为[STIR]中约束边界最小化问题的一阶最优性条件。该算法通过迭代解决增强的信赖域子问题,其特殊对角二次项增强了与边界的距离和梯度方向决定的信赖域形状。这些增强有助于避免直接朝着边界步进,并有效地探索变量空间的整体。为了进一步提高收敛速度,该算法考虑了从边界反射的搜索方向。为了符合理论要求,算法保持迭代点始终可行。对于稠密雅可比矩阵,信赖域子问题通过一种与[JJMore]描述的方法非常相似的准确方法来解决(在 MINPACK 中实现)。与 MINPACK 实现的不同之处在于,雅可比矩阵的奇异值分解在每次迭代中只进行一次,而不是 QR 分解和一系列 Givens 旋转消元。对于大型稀疏雅可比矩阵,使用二维子空间方法解决信赖域子问题,该子空间由梯度的缩放和scipy.sparse.linalg.lsmr
提供的近似 Gauss-Newton 解决方案所构成[STIR],[Byrd]。当不加约束时,该算法与 MINPACK 非常相似,并且通常具有可比较的性能。该算法在无界和有界问题中表现相当强大,因此被选择为默认算法。
方法‘dogbox’在信赖域框架中运作,但考虑到矩形信赖域,而不是传统的椭球体[Voglis]。当前信赖域与初始边界的交集再次是矩形,因此在每次迭代中,通过 Powell 的 dogleg 方法近似地解决了一个受限二次最小化问题[NumOpt]。对于稠密雅可比矩阵,可以通过准确计算 Gauss-Newton 步骤,或者对于大型稀疏雅可比矩阵,可以通过scipy.sparse.linalg.lsmr
进行近似计算。当雅可比矩阵的秩小于变量数时,该算法可能展现出较慢的收敛速度。在变量数较少的有界问题中,该算法通常优于‘trf’。
鲁棒损失函数的实现如 [BA] 描述的那样。其思想是在每次迭代中修改残差向量和雅可比矩阵,使得计算的梯度和 Gauss-Newton Hessian 近似与成本函数的真实梯度和 Hessian 近似匹配。然后算法按正常方式继续,即,鲁棒损失函数被实现为标准最小二乘算法的简单包装。
新版本 0.17.0 中的新增内容。
参考文献
[STIR] (1,2,3)
M. A. Branch, T. F. Coleman 和 Y. Li,《A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems》,《SIAM Journal on Scientific Computing》,第 21 卷,第 1 期,pp. 1-23,1999。
[牛顿拉夫逊法]
William H. Press 等人,《Numerical Recipes. The Art of Scientific Computing. 3rd edition》,第 5.7 节。
[Byrd] (1,2)
R. H. Byrd, R. B. Schnabel 和 G. A. Shultz,《Approximate solution of the trust region problem by minimization over two-dimensional subspaces》,《Math. Programming》,第 40 卷,pp. 247-263,1988。
[Curtis]
A. Curtis, M. J. D. Powell 和 J. Reid,《On the estimation of sparse Jacobian matrices》,《Journal of the Institute of Mathematics and its Applications》,第 13 卷,pp. 117-120,1974。
[JJMore] (1,2,3)
J. J. More,《The Levenberg-Marquardt Algorithm: Implementation and Theory》,《Numerical Analysis》,编者 G. A. Watson,Lecture Notes in Mathematics 630,Springer Verlag,pp. 105-116,1977。
[Voglis]
C. Voglis 和 I. E. Lagaris,《A Rectangular Trust Region Dogleg Approach for Unconstrained and Bound Constrained Nonlinear Optimization》,WSEAS International Conference on Applied Mathematics,希腊科尔富,2004。
[数值优化]
J. Nocedal 和 S. J. Wright,《Numerical optimization, 2nd edition》,第四章。
[BA]
B. Triggs 等人,《Bundle Adjustment - A Modern Synthesis》,《Proceedings of the International Workshop on Vision Algorithms: Theory and Practice》,pp. 298-372,1999。
示例
在这个例子中,我们在独立变量上找到了 Rosenbrock 函数的最小值,没有约束。
>>> import numpy as np
>>> def fun_rosenbrock(x):
... return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
注意,我们只提供残差向量。算法将成本函数构建为残差平方和,这给出了 Rosenbrock 函数。精确的最小值在 x = [1.0, 1.0]
处。
>>> from scipy.optimize import least_squares
>>> x0_rosenbrock = np.array([2, 2])
>>> res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
>>> res_1.x
array([ 1., 1.])
>>> res_1.cost
9.8669242910846867e-30
>>> res_1.optimality
8.8928864934219529e-14
现在我们限制变量,使得先前的解决方案变得不可行。具体来说,我们要求 x[1] >= 1.5
,而 x[0]
保持未约束。为此,我们在 least_squares
中的 bounds 参数中指定了形式为 bounds=([-np.inf, 1.5], np.inf)
。
我们还提供了分析雅可比矩阵:
>>> def jac_rosenbrock(x):
... return np.array([
... [-20 * x[0], 10],
... [-1, 0]])
将所有这些放在一起,我们看到新解决方案位于边界上:
>>> res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
... bounds=([-np.inf, 1.5], np.inf))
>>> res_2.x
array([ 1.22437075, 1.5 ])
>>> res_2.cost
0.025213093946805685
>>> res_2.optimality
1.5885401433157753e-07
现在我们解一个方程组(即,成本函数在最小值处应为零),对一个有 100000 个变量的 Broyden 三对角向量值函数进行约束:
>>> def fun_broyden(x):
... f = (3 - x) * x + 1
... f[1:] -= x[:-1]
... f[:-1] -= 2 * x[1:]
... return f
相应的雅可比矩阵是稀疏的。我们告诉算法通过有限差分来估计它,并提供雅可比矩阵的稀疏结构,以显著加快这个过程。
>>> from scipy.sparse import lil_matrix
>>> def sparsity_broyden(n):
... sparsity = lil_matrix((n, n), dtype=int)
... i = np.arange(n)
... sparsity[i, i] = 1
... i = np.arange(1, n)
... sparsity[i, i - 1] = 1
... i = np.arange(n - 1)
... sparsity[i, i + 1] = 1
... return sparsity
...
>>> n = 100000
>>> x0_broyden = -np.ones(n)
...
>>> res_3 = least_squares(fun_broyden, x0_broyden,
... jac_sparsity=sparsity_broyden(n))
>>> res_3.cost
4.5687069299604613e-23
>>> res_3.optimality
1.1650454296851518e-11
让我们还用鲁棒损失函数解决一个曲线拟合问题,以处理数据中的离群值。定义模型函数为y = a + b * exp(c * t)
,其中 t 是预测变量,y 是观测值,a、b、c 是要估计的参数。
首先,定义生成带有噪声和离群值数据的函数,定义模型参数,并生成数据:
>>> from numpy.random import default_rng
>>> rng = default_rng()
>>> def gen_data(t, a, b, c, noise=0., n_outliers=0, seed=None):
... rng = default_rng(seed)
...
... y = a + b * np.exp(t * c)
...
... error = noise * rng.standard_normal(t.size)
... outliers = rng.integers(0, t.size, n_outliers)
... error[outliers] *= 10
...
... return y + error
...
>>> a = 0.5
>>> b = 2.0
>>> c = -1
>>> t_min = 0
>>> t_max = 10
>>> n_points = 15
...
>>> t_train = np.linspace(t_min, t_max, n_points)
>>> y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)
定义计算残差和参数初始估计的函数。
>>> def fun(x, t, y):
... return x[0] + x[1] * np.exp(x[2] * t) - y
...
>>> x0 = np.array([1.0, 1.0, 0.0])
计算标准最小二乘解:
>>> res_lsq = least_squares(fun, x0, args=(t_train, y_train))
现在使用两个不同的鲁棒损失函数计算两个解。参数f_scale设为 0.1,意味着内点残差不应显著超过 0.1(所用的噪声水平)。
>>> res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,
... args=(t_train, y_train))
>>> res_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,
... args=(t_train, y_train))
最后,绘制所有曲线。我们可以看到通过选择适当的loss,即使在存在强离群值的情况下,我们也可以得到接近最优的估计。但请记住,通常建议首先尝试使用‘soft_l1’或‘huber’损失(如有必要),因为其他两个选项可能会在优化过程中造成困难。
>>> t_test = np.linspace(t_min, t_max, n_points * 10)
>>> y_true = gen_data(t_test, a, b, c)
>>> y_lsq = gen_data(t_test, *res_lsq.x)
>>> y_soft_l1 = gen_data(t_test, *res_soft_l1.x)
>>> y_log = gen_data(t_test, *res_log.x)
...
>>> import matplotlib.pyplot as plt
>>> plt.plot(t_train, y_train, 'o')
>>> plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
>>> plt.plot(t_test, y_lsq, label='linear loss')
>>> plt.plot(t_test, y_soft_l1, label='soft_l1 loss')
>>> plt.plot(t_test, y_log, label='cauchy loss')
>>> plt.xlabel("t")
>>> plt.ylabel("y")
>>> plt.legend()
>>> plt.show()
在下一个示例中,我们展示了如何使用least_squares()
来优化复变量的复值残差函数。考虑以下函数:
>>> def f(z):
... return z - (0.5 + 0.5j)
我们将其封装为一个处理实部和虚部作为独立变量的实变量函数,返回实残差:
>>> def f_wrap(x):
... fx = f(x[0] + 1j*x[1])
... return np.array([fx.real, fx.imag])
因此,我们不再优化原始的 n 个复变量的 m-D 复杂函数,而是优化一个 2n 个实变量的 2m-D 实函数:
>>> from scipy.optimize import least_squares
>>> res_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))
>>> z = res_wrapped.x[0] + res_wrapped.x[1]*1j
>>> z
(0.49999999999925893+0.49999999999925893j)
scipy.optimize.nnls
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.nnls.html#scipy.optimize.nnls
scipy.optimize.nnls(A, b, maxiter=None, *, atol=None)
求解 argmin_x || Ax - b ||_2
使得 x>=0
。
这个问题通常称为非负最小二乘问题,是一个具有凸约束的凸优化问题。当 x
模型的数量只能取得非负值时,通常出现在成分重量、组件成本等方面。
参数:
A(m, n) ndarray
系数数组
b(m,) ndarray, float
右手边向量。
maxiter: int, optional
最大迭代次数,可选。默认值是 3 * n
。
atol: float
在算法中用于评估投影残差 (A.T @ (A x - b)
条目接近零的容差值。增加此值可以放宽解的约束条件。可以选择的典型放宽值为 max(m, n) * np.linalg.norm(a, 1) * np.spacing(1.)
。由于大问题的规范运算变得昂贵,因此此值不设置为默认值,仅在必要时使用。
返回:
xndarray
解向量。
rnormfloat
残差的二范数,|| Ax-b ||_2
。
参见
具有变量界限的线性最小二乘
笔记
该代码基于[2],它是[1]经典算法的改进版本。它利用主动集方法,并解决非负最小二乘问题的 KKT(Karush-Kuhn-Tucker)条件。
参考文献
- [1]
- C. Lawson, R.J. Hanson,“解最小二乘问题”,SIAM,1995,DOI:10.1137/1.9781611971217
- [2]
- Rasmus Bro, Sijmen de Jong,“一种快速的非负约束最小二乘算法”,化学计量学杂志,1997,DOI:10.1002/(SICI)1099-128X(199709/10)11:5<393::AID-CEM483>3.0.CO;2-L
示例
>>> import numpy as np
>>> from scipy.optimize import nnls
...
>>> A = np.array([[1, 0], [1, 0], [0, 1]])
>>> b = np.array([2, 1, 1])
>>> nnls(A, b)
(array([1.5, 1\. ]), 0.7071067811865475)
>>> b = np.array([-1, -1, -1])
>>> nnls(A, b)
(array([0., 0.]), 1.7320508075688772)
scipy.optimize.lsq_linear
scipy.optimize.lsq_linear(A, b, bounds=(-inf, inf), method='trf', tol=1e-10, lsq_solver=None, lsmr_tol=None, max_iter=None, verbose=0, *, lsmr_maxiter=None)
解决具有变量界限的线性最小二乘问题。
给定一个 m × n 的设计矩阵 A 和一个具有 m 个元素的目标向量 b,lsq_linear
解决以下优化问题:
minimize 0.5 * ||A x - b||**2
subject to lb <= x <= ub
该优化问题是凸的,因此找到的最小值(如果迭代收敛)保证是全局的。
参数:
Aarray_like,稀疏矩阵或 LinearOperator,形状为 (m, n)
设计矩阵。可以是 scipy.sparse.linalg.LinearOperator
。
barray_like,形状为 (m,)
目标向量。
bounds2-tuple of array_like 或 Bounds
的实例,可选
参数的上下界。默认情况下没有界限。有两种指定界限的方式:
Bounds
类的实例。- 2-tuple of array_like:元组的每个元素必须是长度等于参数数目的数组,或者是一个标量(在这种情况下,界限被认为对所有参数都是相同的)。使用
np.inf
和适当的符号来禁用所有或某些参数的界限。
method‘trf’ 或 ‘bvls’,可选
执行最小化的方法。
- ‘trf’:适用于线性最小二乘问题的信任区域反射算法。这是一种类似内点的方法,所需的迭代次数与变量数目弱相关。
- ‘bvls’:有界变量最小二乘算法。这是一种活动集方法,需要的迭代次数与变量数目相当。当 A 是稀疏的或 LinearOperator 时无法使用。
默认为 ‘trf’。
tolfloat,可选
容差参数。如果成本函数的相对变化在最后一次迭代中小于 tol,则算法终止。此外,还考虑第一阶优化度量:
method='trf'
如果梯度的均匀范数(考虑到界限的存在)小于 tol,则终止。method='bvls'
如果在 tol 的容差内满足 Karush-Kuhn-Tucker 条件,则终止。
lsq_solver,可选
在迭代过程中解决无界最小二乘问题的方法:
- ‘exact’:使用密集的 QR 或 SVD 分解方法。当 A 是稀疏的或 LinearOperator 时无法使用。
- ‘lsmr’:使用
scipy.sparse.linalg.lsmr
迭代过程,仅需要矩阵-向量乘积评估。不能与method='bvls'
同时使用。
如果为 None(默认值),则根据 A 的类型选择求解器。
lsmr_tolNone、float 或 ‘auto’,可选
耐受参数 ‘atol’ 和 ‘btol’ 用于 scipy.sparse.linalg.lsmr
。如果为 None(默认值),则设置为 1e-2 * tol
。如果是 ‘auto’,则基于当前迭代的最优性调整容差,这可以加速优化过程,但不总是可靠。
max_iterNone 或 int,可选
终止前的最大迭代次数。如果为 None(默认值),则对于 method='trf'
设置为 100,对于 method='bvls'
设置为变量数(不计算 ‘bvls’ 初始化的迭代)。
verbose,可选
算法详细程度:
- 0:静默工作(默认值)。
- 1:显示终止报告。
- 2:显示迭代过程。
lsmr_maxiterNone 或 int,可选
lsmr 最小二乘求解器的最大迭代次数(通过设置 lsq_solver='lsmr'
)。如果为 None(默认值),则使用 lsmr 的默认值 min(m, n)
,其中 m
和 n
分别为 A 的行数和列数。如果 lsq_solver='exact'
,则不起作用。
返回:
OptimizeResult,其以下字段已定义:
xndarray,形状为 (n,)
找到解。
costfloat
解处的成本函数值。
fun数组,形状为 (m,)
解处的残差向量。
optimalityfloat
一阶优化度量。确切含义取决于 method,请参阅 tol 参数的描述。
active_maskint 数组,形状为 (n,)
每个组件显示相应约束是否活跃(即变量是否位于边界):
- 0:无约束被激活。
- -1:下限被激活。
- 1:上限被激活。
对于 trf 方法可能有些随意,因为它生成严格可行迭代序列,并且在容差阈值内确定 active_mask。
unbounded_sol元组
最小二乘求解器返回的无界解元组(使用 lsq_solver 选项设置)。如果 lsq_solver 未设置或设置为 'exact'
,则元组包含形状为 (n,) 的 ndarray,其无界解、残差平方和的 ndarray、A 的秩和 A 的奇异值的 int(请参阅 NumPy 的 linalg.lstsq
获取更多信息)。如果 lsq_solver 设置为 'lsmr'
,则元组包含形状为 (n,) 的 ndarray,其无界解、退出代码的 int、迭代次数的 int 和五个不同规范及 A 的条件数的 float(请参阅 SciPy 的 sparse.linalg.lsmr
获取更多信息)。此输出对于确定最小二乘求解器的收敛性尤为有用,特别是迭代 'lsmr'
求解器。无界最小二乘问题是最小化 0.5 * ||A x - b||**2
。
nitint
迭代次数。如果无约束解是最优解,则为零。
statusint
算法终止的原因:
- -1:算法在最后一次迭代时无法取得进展。
- 0:超过最大迭代次数。
- 1:一阶优化性能度量小于tol。
- 2:成本函数的相对变化小于tol。
- 3:无约束解决方案是最优的。
message字符串
终止原因的口头描述。
success布尔值
如果满足一个收敛标准(status > 0),则为真。
另请参阅
nnls
具有非负约束的线性最小二乘法。
least_squares
具有变量界限的非线性最小二乘法。
注释
该算法首先通过numpy.linalg.lstsq
或scipy.sparse.linalg.lsmr
计算无约束最小二乘解决方案,具体取决于lsq_solver。如果解决方案在界限内,则返回此解决方案作为最优解。
方法 'trf' 运行了适应于线性最小二乘问题的算法描述的改编版[STIR]。迭代基本与非线性最小二乘算法相同,但由于二次函数模型始终准确,因此我们不需要跟踪或修改信任区域的半径。当所选步骤未减少成本函数时,使用线搜索(回溯)作为安全网。详细了解该算法的更多信息,请参阅scipy.optimize.least_squares
。
方法 'bvls' 运行了一个 Python 实现的算法,描述在[BVLS]。该算法维护变量的活动和自由集,在每次迭代中选择一个新变量从活动集移动到自由集,然后在自由变量上解决无约束最小二乘问题。此算法保证最终提供准确的解决方案,但对于具有 n 个变量的问题可能需要多达 n 次迭代。此外,还实施了一种特定初始化过程,确定最初要设置为自由或活动的变量。在实际 BVLS 开始之前需要进行一些迭代,但可以显著减少进一步迭代次数。
参考文献
[STIR]
M. A. Branch, T. F. Coleman 和 Y. Li,《大规模约束最小化问题的子空间、内点和共轭梯度法》,《SIAM 科学计算杂志》,第 21 卷,第 1 号,1-23 页,1999 年。
[BVLS]
P. B. Start 和 R. L. Parker,《有界变量最小二乘法:算法与应用》,《计算统计学》,10,129-141,1995 年。
示例
在这个例子中,解决了一个涉及大稀疏矩阵和变量边界的问题。
>>> import numpy as np
>>> from scipy.sparse import rand
>>> from scipy.optimize import lsq_linear
>>> rng = np.random.default_rng()
...
>>> m = 20000
>>> n = 10000
...
>>> A = rand(m, n, density=1e-4, random_state=rng)
>>> b = rng.standard_normal(m)
...
>>> lb = rng.standard_normal(n)
>>> ub = lb + 1
...
>>> res = lsq_linear(A, b, bounds=(lb, ub), lsmr_tol='auto', verbose=1)
# may vary
The relative change of the cost function is less than `tol`.
Number of iterations 16, initial cost 1.5039e+04, final cost 1.1112e+04,
first-order optimality 4.66e-08.
scipy.optimize.isotonic_regression
scipy.optimize.isotonic_regression(y, *, weights=None, increasing=True)
非参数等温回归。
通过池相邻违反者算法(PAVA)计算出与y长度相同的(不严格)单调递增数组x,参见[1]。更多细节请参见注释部分。
参数:
y(N,) array_like
响应变量。
weights(N,) array_like or None
案例权重。
increasingbool
如果为 True,则拟合单调递增,即等温,回归。如果为 False,则拟合单调递减,即反等温,回归。默认为 True。
返回:
resOptimizeResult
优化结果表示为OptimizeResult
对象。重要属性包括:
-
x
:等温回归解,即与 y 长度相同的递增(或递减)数组,元素范围从 min(y)到 max(y)。 -
weights
:每个块(或池)B 的案例权重总和的数组。 -
blocks
:长度为 B+1 的数组,其中包含每个块(或池)B 的起始位置的索引。第 j 个块由x[blocks[j]:blocks[j+1]]
给出,其中所有值都相同。
注释
给定数据(y)和案例权重(w),等温回归解决了以下优化问题:
[\operatorname{argmin}_{x_i} \sum_i w_i (y_i - x_i)² \quad \text{subject to } x_i \leq x_j \text{ whenever } i \leq j ,.]
对于每个输入值(y_i),它生成一个值(x_i),使得(x)是递增的(但不是严格的),即(x_i \leq x_{i+1})。这是通过 PAVA 完成的。解决方案由池或块组成,即(x)的相邻元素,例如(x_i)和(x_{i+1}),它们都具有相同的值。
最有趣的是,如果将平方损失替换为广泛的 Bregman 函数类,那么解决方案将保持不变,这些函数是均值的唯一一类严格一致的评分函数,参见[2]及其中的参考文献。
根据[1]实现的 PAVA 版本,其计算复杂度为 O(N),其中 N 为输入大小。
参考文献
[1] (1,2)
Busing, F. M. T. A. (2022). 单调回归:简单快速的 O(n) PAVA 实现。《统计软件杂志》,代码片段,102(1),1-25。DOI:10.18637/jss.v102.c01
[2]
Jordan, A.I., Mühlemann, A. & Ziegel, J.F. 表征可识别函数的等温回归问题的最优解。《统计数学研究所通报》74,489-514 (2022)。DOI:10.1007/s10463-021-00808-0
示例
该示例演示了isotonic_regression
确实解决了一个受限制的优化问题。
>>> import numpy as np
>>> from scipy.optimize import isotonic_regression, minimize
>>> y = [1.5, 1.0, 4.0, 6.0, 5.7, 5.0, 7.8, 9.0, 7.5, 9.5, 9.0]
>>> def objective(yhat, y):
... return np.sum((yhat - y)**2)
>>> def constraint(yhat, y):
... # This is for a monotonically increasing regression.
... return np.diff(yhat)
>>> result = minimize(objective, x0=y, args=(y,),
... constraints=[{'type': 'ineq',
... 'fun': lambda x: constraint(x, y)}])
>>> result.x
array([1.25 , 1.25 , 4\. , 5.56666667, 5.56666667,
5.56666667, 7.8 , 8.25 , 8.25 , 9.25 ,
9.25 ])
>>> result = isotonic_regression(y)
>>> result.x
array([1.25 , 1.25 , 4\. , 5.56666667, 5.56666667,
5.56666667, 7.8 , 8.25 , 8.25 , 9.25 ,
9.25 ])
相对于调用minimize
,isotonic_regression
的一个巨大优势在于它更加用户友好,即无需定义目标和约束函数,并且速度快上几个数量级。在普通硬件(2023 年)上,对长度为 1000 的正态分布输入 y 进行优化,最小化器大约需要 4 秒,而isotonic_regression
只需大约 200 微秒。
scipy.optimize.curve_fit
scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=None, bounds=(-inf, inf), method=None, jac=None, *, full_output=False, nan_policy=None, **kwargs)
使用非线性最小二乘拟合函数 f 到数据。
假设 ydata = f(xdata, *params) + eps
。
参数:
f:callable
模型函数,f(x, …)。它必须将独立变量作为第一个参数,将要拟合的参数作为单独的剩余参数。
xdata:array_like
数据测量的自变量。通常应为长度为 M 的序列或形状为 (k,M) 的数组,对于具有 k 个预测变量的函数,如果是类似数组的对象,则每个元素应该可转换为 float。
ydata:array_like
依赖数据,长度为 M 的数组 - 名义上 f(xdata, ...)
。
p0:array_like,可选
参数的初始猜测(长度为 N)。如果为 None,则所有初始值将为 1(如果可以使用内省确定函数的参数数量,否则会引发 ValueError)。
sigma:None 或标量或长度为 M 的序列或 MxM 数组,可选
确定 ydata 中的不确定性。如果定义残差为 r = ydata - f(xdata, *popt)
,那么 sigma 的解释取决于它的维数:
- 一个标量或 1-D sigma 应包含 ydata 中误差的标准偏差值。在这种情况下,优化的函数为
chisq = sum((r / sigma) ** 2)
。- 一个 2-D sigma 应包含 ydata 中误差的协方差矩阵。在这种情况下,优化的函数为
chisq = r.T @ inv(sigma) @ r
。- 新版本为 0.19。
None(默认)等效于填充为 1 的 1-D sigma。
absolute_sigma:bool,可选
如果为 True,则 sigma 以绝对意义使用,并且估计的参数协方差 pcov 反映这些绝对值。
如果为 False(默认),则仅相对大小的 sigma 值有关。返回的参数协方差矩阵 pcov 是通过将 sigma 缩放一个常数因子来计算的。这个常数是通过要求在使用 scaled sigma 时,最优参数 popt 的减少的 chisq 等于单位来设定的。换句话说,sigma 被缩放以匹配拟合后残差的样本方差。默认为 False。数学上,pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)
check_finite:bool,可选
如果为 True,则检查输入数组是否不包含 nans 或 infs,并在包含时引发 ValueError。如果输入数组包含 nans,则将此参数设置为 False 可能会无声地产生荒谬的结果。如果 nan_policy 未明确指定,则默认为 True,否则为 False。
bounds:2 元组的 array_like 或 Bounds
参数的下界和上界。默认无边界。有两种指定边界的方法:
Bounds
类的实例。- 2-tuple 的 array_like:元组的每个元素必须是与参数数量相等的长度的 array 或标量(在这种情况下,边界被视为对所有参数相同)。使用
np.inf
和适当的符号来禁用所有或部分参数的边界。
method {‘lm’, ‘trf’, ‘dogbox’},可选
优化使用的方法。详见 least_squares
以获取更多细节。默认为 ‘lm’ 用于无约束问题和 ‘trf’ 如果提供了 bounds。当观测数量少于变量数量时,方法 ‘lm’ 将无法工作,此时使用 ‘trf’ 或 ‘dogbox’。
新功能在版本 0.17 中引入。
jac 可调用函数、字符串或 None,可选
带有签名 jac(x, ...)
的函数,计算模型函数相对于参数的雅可比矩阵作为密集的 array_like 结构。它将根据提供的 sigma 进行缩放。如果为 None(默认),则将数值地估计雅可比矩阵。可以使用 ‘trf’ 和 ‘dogbox’ 方法的字符串关键字来选择有限差分方案,请参阅 least_squares
。
新功能在版本 0.18 中引入。
full_output 布尔值,可选
如果为 True,此函数将返回额外的信息:infodict、mesg 和 ier。
新功能在版本 1.9 中引入。
nan_policy {‘raise’, ‘omit’, None},可选
定义当输入包含 NaN 时如何处理。可用以下选项(默认为 None):
- ‘raise’:抛出一个错误
- ‘omit’:在计算时忽略 NaN 值
- None:不执行 NaN 的特殊处理(除了 check_finite 执行的内容);当存在 NaN 时的行为取决于实现,并且可能会更改。
请注意,如果显式指定了此值(而不是 None),check_finite 将设置为 False。
新功能在版本 1.11 中引入。
**kwargs
传递给 method='lm'
的 leastsq
或否则传递给 least_squares
的关键字参数。
返回:
popt 数组
优化参数的最佳值,以使 f(xdata, *popt) - ydata
的残差平方和最小化。
pcov 2-D 数组
popt 的估计近似协方差。对角线提供参数估计的方差。要计算参数的一标准差误差,使用 perr = np.sqrt(np.diag(pcov))
。注意 cov 与参数误差估计之间的关系是基于模型函数在最优解周围的线性近似 [1]。当此近似不准确时,cov 可能不提供准确的不确定性测量。
sigma参数如何影响估计协方差取决于absolute_sigma参数,如上所述。
如果解的雅可比矩阵没有完全秩,则‘lm’方法返回一个填满np.inf
的矩阵,而‘trf’和‘dogbox’方法则使用 Moore-Penrose 伪逆来计算协方差矩阵。具有大条件数的协方差矩阵(例如使用numpy.linalg.cond
计算的协方差矩阵)可能表明结果不可靠。
infodictdict(仅在full_output为 True 时返回)
一个带有键的可选输出字典:
nfev
函数调用次数。方法‘trf’和‘dogbox’不对数值雅可比逼近计数函数调用,而‘lm’方法则计算。
fvec
在解决方案处评估的残差值,对于 1-D sigma,这是(f(x, *popt) - ydata)/sigma
。
fjac
一个 QR 因子分解的 R 矩阵的列置换,以列顺序存储。与 ipvt 一起,可以近似估计估计的协方差。‘lm’方法仅提供此信息。
ipvt
长度为 N 的整数数组,定义一个置换矩阵 p,使得 fjacp = qr,其中 r 是对角元素非递增的上三角形矩阵。p 的第 j 列是单位矩阵的 ipvt(j)列。‘lm’方法仅提供此信息。
qtf
向量(转置(q) * fvec)。‘lm’方法仅提供此信息。
新版本 1.9 中的新功能。
mesgstr(仅在full_output为 True 时返回)
一个提供关于解决方案信息的字符串消息。
新版本 1.9 中的新功能。
ierint(仅在full_output为 True 时返回)
一个整数标志。如果等于 1、2、3 或 4,则找到了解决方案。否则,未找到解决方案。在任何情况下,可选输出变量mesg提供更多信息。
新版本 1.9 中的新功能。
Raises:
ValueError
如果ydata或xdata包含 NaN,或者使用不兼容的选项。
RuntimeError
如果最小二乘法最小化失败。
OptimizeWarning
如果无法估计参数的协方差。
另请参阅
least_squares
最小化非线性函数的平方和。
scipy.stats.linregress
为两组测量计算线性最小二乘回归。
注释
用户应确保输入xdata、ydata和f的输出为float64
,否则优化可能返回不正确的结果。
使用method='lm'
,算法通过leastsq
使用 Levenberg-Marquardt 算法。请注意,此算法只能处理无约束问题。
箱约束可以通过‘trf’和‘dogbox’方法处理。有关更多信息,请参阅least_squares
的文档字符串。
参考文献
[1] K. Vugrin 等。非线性置信区间估计技术
地下水流回归:三个案例研究。水资源研究,第 43 卷,W03423,DOI:10.1029/2005WR004804
示例
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.optimize import curve_fit
>>> def func(x, a, b, c):
... return a * np.exp(-b * x) + c
定义要拟合的带有一些噪声的数据:
>>> xdata = np.linspace(0, 4, 50)
>>> y = func(xdata, 2.5, 1.3, 0.5)
>>> rng = np.random.default_rng()
>>> y_noise = 0.2 * rng.normal(size=xdata.size)
>>> ydata = y + y_noise
>>> plt.plot(xdata, ydata, 'b-', label='data')
对函数func的参数 a、b、c 进行拟合:
>>> popt, pcov = curve_fit(func, xdata, ydata)
>>> popt
array([2.56274217, 1.37268521, 0.47427475])
>>> plt.plot(xdata, func(xdata, *popt), 'r-',
... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
优化约束在区域0 <= a <= 3
,0 <= b <= 1
和0 <= c <= 0.5
:
>>> popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5]))
>>> popt
array([2.43736712, 1\. , 0.34463856])
>>> plt.plot(xdata, func(xdata, *popt), 'g--',
... label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))
>>> plt.xlabel('x')
>>> plt.ylabel('y')
>>> plt.legend()
>>> plt.show()
为了可靠的结果,模型func不应该过于参数化;多余的参数可能导致不可靠的协方差矩阵,并且在某些情况下,拟合质量较差。作为对模型是否过于参数化的快速检查,计算协方差矩阵的条件数:
>>> np.linalg.cond(pcov)
34.571092161547405 # may vary
值很小,所以并不引起太多关注。然而,如果我们要向func添加第四个参数d
,其效果与a
相同:
>>> def func(x, a, b, c, d):
... return a * d * np.exp(-b * x) + c # a and d are redundant
>>> popt, pcov = curve_fit(func, xdata, ydata)
>>> np.linalg.cond(pcov)
1.13250718925596e+32 # may vary
这样一个大的值是令人担忧的。协方差矩阵的对角线元素与拟合不确定性相关,提供更多信息:
>>> np.diag(pcov)
array([1.48814742e+29, 3.78596560e-02, 5.39253738e-03, 2.76417220e+28]) # may vary
请注意,第一个和最后一个术语远大于其他元素,表明这些参数的最优值是不明确的,模型中只需要其中一个参数。
scipy.optimize.root_scalar
scipy.optimize.root_scalar(f, args=(), method=None, bracket=None, fprime=None, fprime2=None, x0=None, x1=None, xtol=None, rtol=None, maxiter=None, options=None)
找到标量函数的根。
参数:
f可调用对象
用于查找根的函数。
args元组,可选
传递给目标函数及其导数的额外参数。
method字符串,可选
求解器的类型。应为以下之一
- ‘bisect’ (请见此处)
- ‘brentq’ (请见此处)
- ‘brenth’ (请见此处)
- ‘ridder’ (请见此处)
- ‘toms748’ (请见此处)
- ‘newton’ (请见此处)
- ‘secant’ (请见此处)
- ‘halley’ (请见此处)
bracket: 两个浮点数的序列,可选
围绕根的区间。*f(x, args)在两个端点处具有不同的符号。
x0浮点数,可选
初始猜测。
x1浮点数,可选
第二个猜测。
fprime布尔值或可调用对象,可选
如果fprime是布尔值并且为 True,则假定f返回目标函数的值及其导数。fprime也可以是一个可调用函数,返回f的导数。在这种情况下,它必须接受与f相同的参数。
fprime2布尔值或可调用对象,可选
如果fprime2是布尔值且为 True,则假定f返回目标函数及其一阶和二阶导数的值。fprime2也可以是一个可调用函数,返回f的二阶导数。在这种情况下,它必须接受与f相同的参数。
xtol浮点数,可选
终止的容忍度(绝对)。
rtol浮点数,可选
终止的容忍度(相对)。
maxiter整数,可选
最大迭代次数。
options字典,可选
求解器选项的字典。例如,k
,详见show_options()
。
返回:
solRootResults
以RootResults
对象表示的解。重要属性包括:root
解,converged
表示算法是否成功退出的布尔标志,flag
描述终止原因。详见RootResults
了解其他属性的描述。
另请参见
show_options
求解器接受的额外选项。
根
找到向量函数的根。
注释
本节描述了可以通过‘方法’参数选择的可用求解器。
默认情况下,将使用适合当前情况的最佳方法。如果提供了一个区间,可能会使用其中一种区间方法。如果指定了导数和初始值,可能会选择其中一种基于导数的方法。如果判断没有适用的方法,将会引发异常。
每种方法的参数如下(x=必须, o=可选)。
方法 | f | args | 区间 | x0 | x1 | fprime | fprime2 | xtol | rtol | 最大迭代次数 | 选项 |
---|---|---|---|---|---|---|---|---|---|---|---|
二分法 | x | o | x | o | o | o | o | ||||
布伦特法 | x | o | x | o | o | o | o | ||||
布伦特-史密斯法 | x | o | x | o | o | o | o | ||||
里德法 | x | o | x | o | o | o | o | ||||
TOMS748 | x | o | x | o | o | o | o | ||||
割线法 | x | o | x | o | o | o | o | o | |||
牛顿法 | x | o | x | o | o | o | o | o | |||
哈雷法 | x | o | x | x | x | o | o | o | o |
示例
找到简单三次函数的根
>>> from scipy import optimize
>>> def f(x):
... return (x**3 - 1) # only one real root at x = 1
>>> def fprime(x):
... return 3*x**2
布伦特法
方法以一个区间作为输入
>>> sol = optimize.root_scalar(f, bracket=[0, 3], method='brentq')
>>> sol.root, sol.iterations, sol.function_calls
(1.0, 10, 11)
牛顿法
方法以单个点作为输入,并使用其导数。
>>> sol = optimize.root_scalar(f, x0=0.2, fprime=fprime, method='newton')
>>> sol.root, sol.iterations, sol.function_calls
(1.0, 11, 22)
该函数可以在单次调用中提供值和导数。
>>> def f_p_pp(x):
... return (x**3 - 1), 3*x**2, 6*x
>>> sol = optimize.root_scalar(
... f_p_pp, x0=0.2, fprime=True, method='newton'
... )
>>> sol.root, sol.iterations, sol.function_calls
(1.0, 11, 11)
>>> sol = optimize.root_scalar(
... f_p_pp, x0=0.2, fprime=True, fprime2=True, method='halley'
... )
>>> sol.root, sol.iterations, sol.function_calls
(1.0, 7, 8)
scipy.optimize.brentq
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq
scipy.optimize.brentq(f, a, b, args=(), xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)
使用 Brent 方法在一个包围区间内找到函数的根。
使用经典的 Brent 方法在符号变化的区间[a, b]上找到函数f的根。通常被认为是这里根查找例程中最好的。它是使用反向二次插值的割线法的安全版本。Brent 方法结合了根的定位、区间二分和反向二次插值。有时也被称为 van Wijngaarden-Dekker-Brent 方法。Brent(1973)声称对[a, b]内可计算函数保证收敛。
[Brent1973]提供了该算法的经典描述。另一个描述可以在最近一版的《Numerical Recipes》中找到,包括[PressEtal1992]。第三种描述位于mathworld.wolfram.com/BrentsMethod.html
。通过阅读我们的代码,应该很容易理解该算法。我们的代码与标准表述有些不同:我们选择了不同的外推步骤公式。
参数:
ffunction
Python 函数返回一个数字。函数(f)必须是连续的,并且(f(a))和(f(b))必须有相反的符号。
ascalar
包围区间([a, b])的一个端点。
bscalar
包围区间([a, b])的另一个端点。
xtolnumber, optional
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x是精确的根。该参数必须是正的。对于良好的函数,Brent 方法通常能满足xtol/2
和rtol/2
的上述条件。[Brent1973]
rtolnumber, optional
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x是精确的根。该参数不能小于其默认值4*np.finfo(float).eps
。对于良好的函数,Brent 方法通常能满足xtol/2
和rtol/2
的上述条件。[Brent1973]
maxiterint, optional
如果在maxiter次迭代中未实现收敛,则会引发错误。必须 >= 0。
argstuple, optional
包含函数f的额外参数。f通过apply(f, (x)+args)
调用。
full_outputbool, optional
如果full_output为 False,则返回根。如果full_output为 True,则返回值是(x, r)
,其中x是根,r是一个RootResults
对象。
dispbool, optional
如果为 True,则在算法未收敛时引发 RuntimeError。否则,收敛状态记录在任何RootResults
返回对象中。
返回:
rootfloat
f在a和b之间的根。
rRootResults
(如果full_output = True
)
包含有关收敛情况的信息对象。特别地,如果例程收敛,则r.converged
为 True。
注意事项
f必须连续。f(a)和 f(b)必须具有相反的符号。
相关函数可分为多个类别:
多元局部优化器
fmin
, fmin_powell
, fmin_cg
, fmin_bfgs
, fmin_ncg
非线性最小二乘最小化器
leastsq
受约束的多元优化器
fmin_l_bfgs_b
, fmin_tnc
, fmin_cobyla
全局优化器
basinhopping
, brute
, differential_evolution
本地标量最小化器
fminbound
, brent
, golden
, bracket
N 维根查找
fsolve
1 维根查找
brenth
, ridder
, bisect
, newton
标量固定点查找器
fixed_point
参考资料
[Brent1973] (1,2,3)
Brent, R. P., 无导数最小化算法. 美国新泽西州恩格尔伍德克利夫斯:Prentice-Hall 出版社,1973 年。第 3-4 章。
[PressEtal1992]
Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; 和 Vetterling, W. T. Numerical Recipes in FORTRAN: 科学计算艺术, 第 2 版。英国剑桥:剑桥大学出版社,1992 年。第 9.3 节:“Van Wijngaarden-Dekker-Brent 方法”。
例子
>>> def f(x):
... return (x**2 - 1)
>>> from scipy import optimize
>>> root = optimize.brentq(f, -2, 0)
>>> root
-1.0
>>> root = optimize.brentq(f, 0, 2)
>>> root
1.0
scipy.optimize.brenth
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.brenth.html#scipy.optimize.brenth
scipy.optimize.brenth(f, a, b, args=(), xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)
使用 Brent 方法和双曲线外推法在括号区间中找到函数的根。
一种变体的经典 Brent 例程,用于在参数 a 和 b 之间找到函数 f 的根,其使用双曲线外推法而不是逆二次外推法。Bus&Dekker(1975)保证了该方法的收敛性,并声称此处的函数评估上限是二分法的 4 或 5 倍。f(a)和 f(b)不能具有相同的符号。通常与 brent 例程相当,但没有经过如此深入的测试。这是一种使用双曲线外推法的安全版本的弦截法。此处的版本由 Chuck Harris 编写,并实现了[BusAndDekker1975]的算法 M,其中可以找到进一步的细节(收敛特性、额外的备注等)。
参数:
f函数
返回一个数字的 Python 函数。f 必须连续,并且 f(a)和 f(b)必须具有相反的符号。
a标量
括号区间的一端[a,b]。
b标量
括号区间的另一端[a,b]。
xtol数字,可选
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。该参数必须为正数。与brentq
一样,对于良好的函数,该方法通常会使用xtol/2
和rtol/2
满足上述条件。
rtol数字,可选
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。该参数不能小于其默认值4*np.finfo(float).eps
。与brentq
一样,对于良好的函数,该方法通常会使用xtol/2
和rtol/2
满足上述条件。
maxiter整数,可选
如果在maxiter次迭代中未达到收敛,则会引发错误。必须 >= 0。
args元组,可选
包含函数f的额外参数。通过apply(f, (x)+args)
调用f。
full_output布尔值,可选
如果full_output为 False,则返回根。如果full_output为 True,则返回值为(x, r)
,其中x是根,r是一个RootResults
对象。
disp布尔值,可选
如果为 True,则如果算法未收敛,则引发 RuntimeError。否则,收敛状态记录在任何RootResults
返回对象中。
返回:
root浮点数
f在a和b之间的根。
rRootResults
(present if full_output = True
)
包含收敛信息的对象。特别地,如果程序收敛,则 r.converged
为 True。
参见
fmin
, fmin_powell
, fmin_cg
, fmin_bfgs
, fmin_ncg
多元局部优化器
leastsq
非线性最小二乘优化器
fmin_l_bfgs_b
, fmin_tnc
, fmin_cobyla
有约束的多元优化器
basinhopping
, differential_evolution
, brute
全局优化器
fminbound
, brent
, golden
, bracket
局部标量最小化器
fsolve
N 元根查找
brentq
, brenth
, ridder
, bisect
, newton
一维根查找
fixed_point
标量固定点查找器
参考文献
[BusAndDekker1975]
Bus, J. C. P., Dekker, T. J., “Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function”, ACM Transactions on Mathematical Software, Vol. 1, Issue 4, Dec. 1975, pp. 330-345. Section 3: “Algorithm M”. DOI:10.1145/355656.355659
示例
>>> def f(x):
... return (x**2 - 1)
>>> from scipy import optimize
>>> root = optimize.brenth(f, -2, 0)
>>> root
-1.0
>>> root = optimize.brenth(f, 0, 2)
>>> root
1.0
scipy.optimize.ridder
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.ridder.html#scipy.optimize.ridder
scipy.optimize.ridder(f, a, b, args=(), xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)
使用 Ridder 方法在区间内查找函数的根。
参数:
f函数
返回一个数字的 Python 函数。f 必须连续,并且 f(a) 和 f(b) 必须有相反的符号。
a标量
区间 [a,b] 的一端。
b标量
区间 [a,b] 的另一端。
xtol数字,可选
计算的根 x0
将满足 np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中 x
是精确的根。参数必须为正。
rtol数字,可选
计算的根 x0
将满足 np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中 x
是精确的根。参数不能小于其默认值 4*np.finfo(float).eps
。
maxiter整数,可选
如果在 maxiter 次迭代中未实现收敛,则会引发错误。必须 >= 0。
args元组,可选
包含用于函数 f 的额外参数。通过 apply(f, (x)+args)
调用 f。
full_output布尔值,可选
如果 full_output 为 False,则返回根。如果 full_output 为 True,则返回 (x, r)
,其中 x 是根,r 是一个 RootResults
对象。
disp布尔值,可选
如果为 True,则在算法未收敛时引发 RuntimeError。否则,收敛状态记录在任何 RootResults
返回对象中。
返回:
root浮点数
f 在 a 和 b 之间的根。
rRootResults
(如果 full_output = True
)
包含有关收敛信息的对象。特别是,如果例程收敛,则 r.converged
为 True。
另请参阅
brentq
,brenth
,bisect
,newton
1-D 根查找
fixed_point
标量固定点查找器
注意
使用 [Ridders1979] 方法在函数 f 的参数 a 和 b 之间找到根。Ridders 方法比二分法更快,但通常不如 Brent 方法快。[Ridders1979] 提供了算法的经典描述和源。在任何最新版本的《数值方法》中也可以找到描述。
此处使用的例行程序略有偏离标准演示,以更加谨慎地处理容差。
参考文献
[Ridders1979] (1,2)
Ridders, C. F. J. “A New Algorithm for Computing a Single Root of a Real Continuous Function.” IEEE Trans. Circuits Systems 26, 979-980, 1979.
示例
>>> def f(x):
... return (x**2 - 1)
>>> from scipy import optimize
>>> root = optimize.ridder(f, 0, 2)
>>> root
1.0
>>> root = optimize.ridder(f, -2, 0)
>>> root
-1.0
scipy.optimize.bisect
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.bisect.html#scipy.optimize.bisect
scipy.optimize.bisect(f, a, b, args=(), xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)
使用二分法在区间内找到函数的根。
基本的二分法例程,用于在参数a和b之间找到函数f的根。f(a)和f(b)不能有相同的符号。缓慢但可靠。
参数:
f函数
返回一个数的 Python 函数。f必须是连续的,且f(a)和f(b)必须有相反的符号。
a标量
一个括号间隔的端点[a,b]。
b标量
括号间隔的另一端[a,b]。
xtol数值,可选
计算得到的根x0
满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。该参数必须为正。
rtol数值,可选
计算得到的根x0
满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。该参数不能小于其默认值4*np.finfo(float).eps
。
maxiter整数,可选
如果在maxiter次迭代中未实现收敛,则引发错误。必须>= 0。
args元组,可选
包含传递给函数f的额外参数。f由apply(f, (x)+args)
调用。
full_output布尔型,可选
如果full_output为 False,则返回根。如果full_output为 True,则返回值为(x, r)
,其中 x 为根,r 为RootResults
对象。
disp布尔型,可选
如果为 True,则在算法未收敛时引发 RuntimeError。否则,收敛状态记录在RootResults
返回对象中。
返回:
root浮点数
f在a和b之间的根。
rRootResults
(如果full_output = True
)
包含有关收敛性的信息的对象。特别地,如果程序收敛,则r.converged
为 True。
另请参阅
标量的固定点查找器
n 维根查找
示例
>>> def f(x):
... return (x**2 - 1)
>>> from scipy import optimize
>>> root = optimize.bisect(f, 0, 2)
>>> root
1.0
>>> root = optimize.bisect(f, -2, 0)
>>> root
-1.0
scipy.optimize.newton
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.newton.html#scipy.optimize.newton
scipy.optimize.newton(func, x0, fprime=None, args=(), tol=1.48e-08, maxiter=50, fprime2=None, x1=None, rtol=0.0, full_output=False, disp=True)
使用牛顿-拉弗森(或割线或哈雷)方法找到实数或复数函数的根。
找到标量值函数func的根,给定附近的标量起始点x0。如果函数func的导数fprime被提供,则使用牛顿-拉弗森方法,否则使用割线法。如果函数func的二阶导数fprime2也被提供,则使用哈雷方法。
如果x0是一个具有多个项的序列,newton
将返回一个数组:从x0中的每个(标量)起始点的函数的根。在这种情况下,func必须被矢量化以返回与其第一个参数相同形状的序列或数组。如果给定fprime(fprime2),则其返回值也必须具有相同的形状:其每个元素是函数func相对于其唯一变量在其第一个参数的每个元素处求值的第一(第二)导数。
newton
用于查找单变量标量函数的根。对于涉及多个变量的问题,请参阅root
。
参数:
funccallable
所需的根函数。它必须是形式为f(x,a,b,c...)
的单变量函数,其中a,b,c...
是可以在args参数中传递的额外参数。
x0float, sequence, or ndarray
应该接近实际根的初始估计。如果不是标量,则func必须被矢量化,并且返回与其第一个参数相同形状的序列或数组。
fprimecallable, optional
当函数的导数可用且方便时。如果为 None(默认),则使用割线法。
argstuple, optional
用于函数调用的额外参数。
tolfloat, optional
根值的允许误差。如果func是复数值的,建议使用较大的tol,因为x的实部和虚部都会影响到|x - x0|
。
maxiterint, optional
最大迭代次数。
fprime2callable, optional
当函数的二阶导数可用且方便时。如果为 None(默认),则使用正常的牛顿-拉弗森或割线法。如果不为 None,则使用哈雷法。
x1float, optional
另一个估计的根,应该接近实际根。如果未提供fprime,则使用。
rtolfloat, optional
终止的容差(相对值)。
full_outputbool, optional
如果full_output为 False(默认),则返回根。如果为 True 并且x0为标量,则返回值为(x, r)
,其中x
是根,r
是RootResults
对象。如果为 True 并且x0为非标量,则返回值为(x, converged, zero_der)
(详见返回部分)。
disp布尔值,可选
如果为 True,则在算法未收敛时引发 RuntimeError,错误消息包含迭代次数和当前函数值。否则,收敛状态记录在RootResults
返回对象中。如果x0不是标量,则忽略。注意:这与显示无关,但为了向后兼容性,不能重命名disp
关键字。
返回:
根浮点数、序列或 ndarray
估计的函数为零的位置。
rRootResults
,可选
如果full_output=True
且x0为标量。包含有关收敛性的信息的对象。特别地,如果例程收敛,则r.converged
为 True。
converged布尔值的 ndarray,可选
如果full_output=True
并且x0为非标量。对于向量函数,指示哪些元素成功收敛。
zero_der布尔值的 ndarray,可选
如果full_output=True
并且x0不是标量。对于向量函数,指示哪些元素具有零导数。
另请参见
root_scalar
标量函数的根求解器接口
根
多输入多输出函数的根求解器接口
注意
牛顿-拉弗森方法的收敛速度是二次的,海莉方法是三次的,割线法是次二次的。这意味着如果函数表现良好,第 n 次迭代后估计根的实际误差大约是第(n-1)步后的平方(海莉方法为立方)。然而,此处使用的停止准则是步长,并不能保证找到根。因此,应验证结果。更安全的算法是 brentq、brenth、ridder 和 bisect,但它们都要求在函数变号的区间中首先找到根。在找到这样的区间后,建议使用 brentq 算法进行一维问题的通用解决。
当使用数组进行newton
时,最适合以下类型的问题:
-
初始猜测值x0相对于根的距离几乎相同。
-
部分或全部的额外参数,args,也是数组,以便可以一起解决一类相似的问题。
-
初始猜测值 x0 的大小大于 O(100) 元素。否则,一个简单的循环可能比向量表现得更好。
示例
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import optimize
>>> def f(x):
... return (x**3 - 1) # only one real root at x = 1
如果只提供了 fprime
,使用割线法:
>>> root = optimize.newton(f, 1.5)
>>> root
1.0000000000000016
>>> root = optimize.newton(f, 1.5, fprime2=lambda x: 6 * x)
>>> root
1.0000000000000016
只有提供了 fprime
,使用牛顿-拉夫逊法:
>>> root = optimize.newton(f, 1.5, fprime=lambda x: 3 * x**2)
>>> root
1.0
如果提供了 fprime2
和 fprime
,使用 Halley 方法:
>>> root = optimize.newton(f, 1.5, fprime=lambda x: 3 * x**2,
... fprime2=lambda x: 6 * x)
>>> root
1.0
当我们想要为一组相关的起始值和/或函数参数找到根时,我们可以将这些作为输入数组提供:
>>> f = lambda x, a: x**3 - a
>>> fder = lambda x, a: 3 * x**2
>>> rng = np.random.default_rng()
>>> x = rng.standard_normal(100)
>>> a = np.arange(-50, 50)
>>> vec_res = optimize.newton(f, x, fprime=fder, args=(a, ), maxiter=200)
上述操作相当于在 for 循环中分别解决每个 (x, a)
值,只是速度更快:
>>> loop_res = [optimize.newton(f, x0, fprime=fder, args=(a0,),
... maxiter=200)
... for x0, a0 in zip(x, a)]
>>> np.allclose(vec_res, loop_res)
True
绘制找到的所有 a
值的结果:
>>> analytical_result = np.sign(a) * np.abs(a)**(1/3)
>>> fig, ax = plt.subplots()
>>> ax.plot(a, analytical_result, 'o')
>>> ax.plot(a, vec_res, '.')
>>> ax.set_xlabel('$a$')
>>> ax.set_ylabel('$x$ where $f(x, a)=0$')
>>> plt.show()
scipy.optimize.toms748
scipy.optimize.toms748(f, a, b, args=(), k=1, xtol=2e-12, rtol=8.881784197001252e-16, maxiter=100, full_output=False, disp=True)
使用 TOMS 算法 748 方法找到根。
实现 Alefeld,Potro 和 Shi 的 Algorithm 748 方法,在区间[a , b]上找到函数f的根,其中f(a)和f(b)必须有相反的符号。
它使用反立方插值和“牛顿二次”步骤的混合。[APS1995]。
参数:
f函数
返回标量的 Python 函数。函数(f)必须连续,并且(f(a))和(f(b))具有相反的符号。
a标量,
搜索区间的下界
b标量,
搜索区间的上界
args元组,可选
包含用于函数f的额外参数的对象。f通过f(x, *args)
调用。
k整数,可选
每次迭代执行的牛顿二次步骤数。k>=1
。
xtol标量,可选
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。该参数必须为正数。
rtol标量,可选
计算得到的根x0
将满足np.allclose(x, x0, atol=xtol, rtol=rtol)
,其中x
是精确的根。
maxiter整数,可选
如果在maxiter次迭代中未收敛,将引发错误。必须大于或等于 0。
full_output布尔值,可选
如果full_output为 False,则返回根。如果full_output为 True,则返回值为(x, r)
,其中x为根,r是RootResults
对象。
disp布尔值,可选
如果为 True,则在算法未收敛时引发运行时错误。否则,收敛状态记录在RootResults
返回对象中。
返回:
root浮点数
f的近似根
rRootResults
(如果full_output = True
时存在)
包含有关收敛性的信息的对象。特别地,如果例程收敛,则r.converged
为 True。
另见
brentq
,brenth
,ridder
,bisect
,newton
在 N 维空间中找到根。
注意
f必须是连续的。算法 748 在k=2
时渐近地是已知寻找四次连续可微函数根最有效的算法。与 Brent 算法相比,在最后一步可能仅减少包围区间的长度,算法 748 在每次迭代中以与找到根的渐近效率相同的方式减小它。
为了便于表述效率指标,假设f具有 4 个连续导数。对于k=1
,收敛阶至少为 2.7,每次迭代约有渐近 2 次函数评估,效率指数约为 1.65。对于k=2
,阶数约为 4.6,每次迭代渐近 3 次函数评估,效率指数为 1.66。对于更高的k值,效率指数接近于(3k-2)
的 k 次根,因此k=1
或k=2
通常是合适的选择。
参考文献
[APS1995]
Alefeld, G. E. and Potra, F. A. and Shi, Yixun,Algorithm 748: Enclosing Zeros of Continuous Functions,ACM Trans. Math. Softw. Volume 221(1995) doi = {10.1145/210089.210111}
示例
>>> def f(x):
... return (x**3 - 1) # only one real root at x = 1
>>> from scipy import optimize
>>> root, results = optimize.toms748(f, 0, 2, full_output=True)
>>> root
1.0
>>> results
converged: True
flag: converged
function_calls: 11
iterations: 5
root: 1.0
method: toms748
scipy.optimize.RootResults
class scipy.optimize.RootResults(root, iterations, function_calls, flag, method)
代表根查找结果。
属性:
rootfloat
估计的根位置。
iterationsint
寻找根所需的迭代次数。
function_callsint
调用函数的次数。
convergedbool
如果例程收敛,则为真。
flagstr
终止原因的描述。
methodstr
使用的根查找方法。
方法
__getitem__ |
x.getitem(y) <==> x[y] |
---|---|
__len__ (/) |
返回 len(self)。 |
clear () |
|
copy () |
|
fromkeys (iterable[, value]) |
使用来自 iterable 的键创建一个新的字典,并将值设置为 value。 |
get (key[, default]) |
如果键在字典中,则返回键的值,否则返回默认值。 |
items () |
|
keys () |
|
pop (key[, default]) |
如果未找到键,则返回给定的默认值,否则引发 KeyError |
popitem (/) |
移除并返回一个(键,值)对作为一个 2 元组。 |
setdefault (key[, default]) |
如果键不在字典中,则插入具有默认值的键。 |
update ([E, ]**F) |
如果 E 存在且具有.keys()方法,则执行以下操作:对于 k 在 E 中:D[k] = E[k] 如果 E 存在但没有.keys()方法,则对于 k,v 在 E 中:D[k] = v 无论哪种情况,接下来进行以下操作:对于 k 在 F 中:D[k] = F[k] |
values () |
scipy.optimize.fixed_point
scipy.optimize.fixed_point(func, x0, args=(), xtol=1e-08, maxiter=500, method='del2')
找到函数的不动点。
给定一个或多个变量的函数和一个起始点,找到函数的不动点,即func(x0) == x0
。
参数:
func函数
评估函数。
x0array_like
函数的不动点。
args元组,可选
func的额外参数。
xtolfloat,可选
收敛容差,默认为 1e-08。
maxiterint,可选
最大迭代次数,默认为 500。
method,可选
找到函数的不动点的方法,默认为“del2”,使用带有 Aitken 的Del²
收敛加速的 Steffensen 方法[1]。"iteration"方法仅迭代函数,直到检测到收敛,而不尝试加速收敛。
参考文献
[1]
Burden, Faires, “Numerical Analysis”, 5th edition, pg. 80
示例
>>> import numpy as np
>>> from scipy import optimize
>>> def func(x, c1, c2):
... return np.sqrt(c1/(x+c2))
>>> c1 = np.array([10,12.])
>>> c2 = np.array([3, 5.])
>>> optimize.fixed_point(func, [1.2, 1.3], args=(c1,c2))
array([ 1.4920333 , 1.37228132])
scipy.optimize.root
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.root.html#scipy.optimize.root
scipy.optimize.root(fun, x0, args=(), method='hybr', jac=None, tol=None, callback=None, options=None)
找到向量函数的根。
参数:
fun可调用对象
用于找到根的向量函数。
x0ndarray
初始猜测。
args元组,可选
传递给目标函数及其雅可比矩阵的额外参数。
methodstr,可选
求解器类型。应为以下之一
- ‘hybr’ (详见此处)
- ‘lm’ (详见此处)
- ‘broyden1’ (详见此处)
- ‘broyden2’ (详见此处)
- ‘anderson’ (详见此处)
- ‘linearmixing’ (详见此处)
- ‘diagbroyden’ (详见此处)
- ‘excitingmixing’ (详见此处)
- ‘krylov’ (详见此处)
- ‘df-sane’ (详见此处)
jacbool 或者可调用对象,可选
如果jac是布尔值且为 True,则fun假定返回雅可比矩阵的值以及目标函数的值。如果为 False,则数值上估计雅可比矩阵。jac也可以是返回fun的雅可比矩阵的可调用对象。在这种情况下,它必须接受与fun相同的参数。
tolfloat,可选
终止容差。要进行详细控制,请使用特定于求解器的选项。
callback函数,可选
可选回调函数。在每次迭代中调用为callback(x, f)
,其中x是当前解,f是相应的残差。适用于所有方法,但不包括‘hybr’和‘lm’。
options字典,可选
求解器选项的字典。例如,xtol或maxiter,详见show_options()
获取详细信息。
返回值:
solOptimizeResult
表示为OptimizeResult
对象的解。重要属性包括:x
解数组,success
表示算法是否成功退出的布尔标志和 message
描述终止原因。详见OptimizeResult
以获取其他属性的描述。
另见
show_options
接受求解器的附加选项
注意
本节描述了可以通过‘method’参数选择的可用求解器。默认方法为hybr。
方法hybr使用的是 MINPACK 中实现的 Powell 混合方法的修改[1]。
方法lm使用修改的 Levenberg-Marquardt 算法解决非线性方程组的最小二乘意义上的问题,其实现在 MINPACK 中。[1]
方法df-sane是一种无导数光谱方法。[3]
方法broyden1、broyden2、anderson、linearmixing、diagbroyden、excitingmixing、krylov都是不精确的牛顿方法,带有回溯或全线搜索[2]。每种方法对应特定的雅可比近似。
-
方法broyden1使用布罗伊登第一雅可比近似,被称为布罗伊登的良好方法。
-
方法broyden2使用布罗伊登第二雅可比近似,被称为布罗伊登的不良方法。
-
方法anderson使用(扩展的)安德森混合。
-
方法Krylov使用克里洛夫逆雅可比近似。适用于大规模问题。
-
方法diagbroyden使用对角布罗伊登雅可比近似。
-
方法linearmixing使用标量雅可比近似。
-
方法excitingmixing使用调整的对角雅可比近似。
警告
为方法diagbroyden、linearmixing和excitingmixing实现的算法可能对特定问题有用,但其是否有效可能极大地依赖于问题本身。
0.11.0 版中的新功能。
参考文献
[1] (1,2)
More, Jorge J., Burton S. Garbow, and Kenneth E. Hillstrom. 1980. MINPACK-1 用户指南。
[2]
C. T. Kelley. 1995. 线性和非线性方程的迭代方法。工业和应用数学学会。<archive.siam.org/books/kelley/fr16/
>
[3]
- La Cruz, J.M. Martinez, M. Raydan. Math. Comp. 75, 1429 (2006).
示例
下列函数定义了非线性方程组及其雅可比矩阵。
>>> import numpy as np
>>> def fun(x):
... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
... 0.5 * (x[1] - x[0])**3 + x[1]]
>>> def jac(x):
... return np.array([[1 + 1.5 * (x[0] - x[1])**2,
... -1.5 * (x[0] - x[1])**2],
... [-1.5 * (x[1] - x[0])**2,
... 1 + 1.5 * (x[1] - x[0])**2]])
可以通过以下方式获得解决方案。
>>> from scipy import optimize
>>> sol = optimize.root(fun, [0, 0], jac=jac, method='hybr')
>>> sol.x
array([ 0.8411639, 0.1588361])
大问题
假设我们需要在正方形([0,1]\times[0,1])上解决以下积分微分方程:
[\nabla² P = 10 \left(\int_0¹\int_0¹\cosh(P),dx,dy\right)²]
满足(P(x,1) = 1)且在正方形边界上(P=0)的其他地方。
可以通过solver='krylov'
找到解决方案:
>>> from scipy import optimize
>>> # parameters
>>> nx, ny = 75, 75
>>> hx, hy = 1./(nx-1), 1./(ny-1)
>>> P_left, P_right = 0, 0
>>> P_top, P_bottom = 1, 0
>>> def residual(P):
... d2x = np.zeros_like(P)
... d2y = np.zeros_like(P)
...
... d2x[1:-1] = (P[2:] - 2*P[1:-1] + P[:-2]) / hx/hx
... d2x[0] = (P[1] - 2*P[0] + P_left)/hx/hx
... d2x[-1] = (P_right - 2*P[-1] + P[-2])/hx/hx
...
... d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
... d2y[:,0] = (P[:,1] - 2*P[:,0] + P_bottom)/hy/hy
... d2y[:,-1] = (P_top - 2*P[:,-1] + P[:,-2])/hy/hy
...
... return d2x + d2y - 10*np.cosh(P).mean()**2
>>> guess = np.zeros((nx, ny), float)
>>> sol = optimize.root(residual, guess, method='krylov')
>>> print('Residual: %g' % abs(residual(sol.x)).max())
Residual: 5.7972e-06 # may vary
>>> import matplotlib.pyplot as plt
>>> x, y = np.mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
>>> plt.pcolormesh(x, y, sol.x, shading='gouraud')
>>> plt.colorbar()
>>> plt.show()
scipy.optimize.milp
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.optimize.milp.html#scipy.optimize.milp
scipy.optimize.milp(c, *, integrality=None, bounds=None, constraints=None, options=None)
混合整数线性规划
解决以下形式的问题:
[\begin{split}\min_x \ & c^T x \ \mbox{使得} \ & b_l \leq A x \leq b_u,\ & l \leq x \leq u, \ & x_i \in \mathbb{Z}, i \in X_i\end{split}]
其中 (x) 是决策变量向量;(c), (b_l), (b_u), (l), 和 (u) 是向量;(A) 是矩阵,(X_i) 是必须是整数的决策变量索引集合。(在此上下文中,只能取整数值的变量称为“整数”;它具有“整数性”约束。)
或者,这样说:
最小化:
c @ x
使得:
b_l <= A @ x <= b_u
l <= x <= u
Specified elements of x must be integers
默认情况下,l = 0
并且 u = np.inf
,除非使用 bounds
进行指定。
参数:
c1D 密集数组
要最小化的线性目标函数的系数。在问题解决之前,c 被转换为双精度数组。
integrality 1D 密集数组,可选
指示每个决策变量的整数约束类型。
0
: 连续变量;无整数约束。
1
: 整数变量;决策变量必须是整数且在边界内。
2
: 半连续变量;决策变量必须在边界内或者取值为 0
。
3
: 半整数变量;决策变量必须是整数且在边界内,或者取值为 0
。
默认情况下,所有变量均为连续变量。整数性在问题解决之前被转换为整数数组。
bounds scipy.optimize.Bounds,可选
决策变量的边界。在问题解决之前,下限和上限被转换为双精度数组。Bounds
对象的 keep_feasible
参数被忽略。如果未指定,则所有决策变量都受到非负约束。
constraints 一系列 scipy.optimize.LinearConstraint,可选
优化问题的线性约束。参数可以是以下之一:
-
单个
LinearConstraint
对象 -
可以转换为单个元组,作为
LinearConstraint
对象的参数LinearConstraint(*constraints)
-
由类型为 1. 和 2. 的对象组成的序列。
在解决问题之前,所有值都转换为双精度,并且约束系数的矩阵转换为scipy.sparse.csc_array
的实例。LinearConstraint
对象的keep_feasible
参数被忽略。
optionsdict,可选
求解器选项的字典。以下键被识别。
dispbool(默认值:False
)
如果要在优化期间将优化状态的指示器打印到控制台,则设置为True
。
node_limitint,可选
解决前停止的最大节点数(线性程序松弛)。默认情况下没有最大节点数。
presolvebool(默认值:True
)
Presolve 尝试在将问题发送给主求解器之前识别微不足道的不可行性,识别微不足道的无界性并简化问题。
time_limitfloat,可选
解决问题的最大秒数。默认情况下没有时间限制。
mip_rel_gapfloat,可选
MIP 求解器的终止准则:当主目标值与对偶目标界限之间的差距,按主目标值缩放,<= mip_rel_gap 时,求解器将终止。
返回:
resOptimizeResult
scipy.optimize.OptimizeResult
的实例。对象保证具有以下属性。
statusint
表示算法退出状态的整数。
0
:找到最优解。
1
:达到迭代或时间限制。
2
:问题不可行。
3
:问题无界。
4
:其他;请参阅详细信息。
successbool
当找到最优解时为True
,否则为False
。
messagestr
算法的退出状态的字符串描述符。
还将存在以下属性,但根据解决方案状态,值可能为None
。
xndarray
决策变量的值,这些值最小化了满足约束条件的目标函数。
funfloat
目标函数c @ x
的最优值。
mip_node_countint
MILP 求解器解决的子问题或“节点”的数量。
mip_dual_boundfloat
MILP 求解器对最优解的下界的最终估计。
mip_gapfloat
主目标值与对偶目标界限之间的差距,按主目标值缩放。
注意事项
milp
是 HiGHS 线性优化软件的包装器[1]。该算法是确定性的,并且通常在存在时找到中度挑战的混合整数线性规划的全局最优解。
参考文献
[1]
Huangfu, Q., Galabova, I., Feldmeier, M., 和 Hall, J. A. J. “HiGHS - 高性能线性优化软件。” highs.dev/
[2]
Huangfu, Q. 和 Hall, J. A. J. “并行化双修正单纯形法。” 数学规划计算, 10 (1), 119-142, 2018. DOI: 10.1007/s12532-017-0130-5
示例
考虑en.wikipedia.org/wiki/Integer_programming#Example
中表达为两个变量最大化问题。由于milp
要求将问题表达为最小化问题,决策变量的目标函数系数为:
>>> import numpy as np
>>> c = -np.array([0, 1])
注意负号:我们通过最小化目标函数的负数来最大化原始目标函数。
我们将约束的系数收集到数组中,例如:
>>> A = np.array([[-1, 1], [3, 2], [2, 3]])
>>> b_u = np.array([1, 12, 12])
>>> b_l = np.full_like(b_u, -np.inf)
由于这些约束没有下限,我们定义了一个变量b_l
,其中包含代表负无穷大的值。这对于scipy.optimize.linprog
的用户可能不熟悉,后者仅接受形式为A_ub @ x <= b_u
的“小于”(或“上界”)不等式约束。通过接受约束b_l <= A_ub @ x <= b_u
的b_l
和b_u
,milp
能够简洁地指定“大于”不等式约束、“小于”不等式约束和等式约束。
将这些数组收集到一个单一的LinearConstraint
对象中,如下:
>>> from scipy.optimize import LinearConstraint
>>> constraints = LinearConstraint(A, b_l, b_u)
决策变量的非负限制默认受到强制执行,因此我们无需为bounds提供参数。
最后,问题规定决策变量必须是整数:
>>> integrality = np.ones_like(c)
我们解决问题如下:
>>> from scipy.optimize import milp
>>> res = milp(c=c, constraints=constraints, integrality=integrality)
>>> res.x
[1.0, 2.0]
请注意,如果我们解决了放松的问题(没有整数约束):
>>> res = milp(c=c, constraints=constraints) # OR:
>>> # from scipy.optimize import linprog; res = linprog(c, A, b_u)
>>> res.x
[1.8, 2.8]
如果我们通过四舍五入到最接近的整数来解决问题,我们将得不到正确的解决方案。
其他示例见于 tutorial 教程。
scipy.optimize.linprog
scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='highs', callback=None, options=None, x0=None, integrality=None)
线性规划:最小化线性目标函数,受线性等式和不等式约束限制。
线性规划解决如下形式的问题:
[\begin{split}\min_x \ & c^T x \ \mbox{使得} \ & A_{ub} x \leq b_{ub},\ & A_{eq} x = b_{eq},\ & l \leq x \leq u ,\end{split}]
其中 (x) 是决策变量向量;(c), (b_{ub}), (b_{eq}), (l), 和 (u) 是向量;(A_{ub}) 和 (A_{eq}) 是矩阵。
或者说:
最小化
c @ x
使得
A_ub @ x <= b_ub A_eq @ x == b_eq lb <= x <= ub
注意,默认情况下 lb = 0
和 ub = None
。可以使用 bounds
指定其他边界。
参数:
c 1-D 数组
要最小化的线性目标函数的系数。
A_ub 2-D 数组,可选
不等约束矩阵。A_ub
的每一行指定 x
的线性不等式约束的系数。
b_ub 1-D 数组,可选
不等约束向量。每个元素表示对应的A_ub @ x
的上限。
A_eq 2-D 数组,可选
等式约束矩阵。A_eq
的每一行指定 x
的线性等式约束的系数。
b_eq 1-D 数组,可选
等式约束向量。A_eq @ x
的每个元素必须等于 b_eq
的对应元素。
bounds 序列,可选
对于 x
中每个元素的 (min, max)
对序列,定义决策变量的最小和最大值。如果提供单个元组 (min, max)
,则 min
和 max
将作为所有决策变量的边界。使用 None
表示无边界。例如,默认边界 (0, None)
表示所有决策变量非负,而对 (None, None)
表示无任何边界,即所有变量可以是任意实数。
method 字符串,可选
用于解决标准形式问题的算法。支持 ‘highs’(默认),‘highs-ds’,‘highs-ipm’,‘interior-point’(遗留),‘revised simplex’(遗留)和 ‘simplex’(遗留)。遗留方法已弃用,将在 SciPy 1.11.0 中移除。
callback 可调用对象,可选
如果提供了回调函数,则算法的每次迭代至少调用一次。回调函数必须接受单一的scipy.optimize.OptimizeResult
,包含以下字段:
x1-D 数组
当前解向量。
funfloat
目标函数 c @ x
的当前值。
successbool
True
when the algorithm has completed successfully.
slack1-D 数组
松弛变量的(名义上的正)值,b_ub - A_ub @ x
。
con1-D 数组
等式约束的(名义上的零)残差,b_eq - A_eq @ x
。
phaseint
正在执行的算法阶段。
statusint
表示算法状态的整数。
0
: 优化正常进行。
1
: 达到迭代限制。
2
: Problem appears to be infeasible.
3
: 问题似乎无界。
4
: 遇到数值困难。
nitint
当前迭代次数。
messagestr
描述算法状态的字符串。
目前 HiGHS 方法不支持回调函数。
optionsdict, optional
求解器选项的字典。所有方法都接受以下选项:
maxiterint
执行的最大迭代次数。默认值:请参阅特定方法的文档。
dispbool
设置为 True
打印收敛消息。默认值:False
。
presolvebool
设置为 False
禁用自动预处理。默认值:True
。
除 HiGHS 求解器外,所有方法都接受:
tolfloat
决定残差何时“足够接近”零以被视为精确零的公差。
autoscalebool
设置为 True
自动执行均衡化。如果约束中的数值在数量级上相隔甚远,则考虑使用此选项。默认值:False
。
rrbool
设置为 False
禁用自动冗余移除。默认值:True
。
rr_methodstring
在预处理后从等式约束矩阵中识别和删除多余行的方法。对于输入稠密的问题,可用的冗余移除方法有:
“SVD”:
反复对矩阵执行奇异值分解,基于左奇异向量中的非零元素检测冗余行,对于接近完全秩的矩阵可能很快。
“pivot”:
使用[5]中介绍的算法来识别多余的行。
“ID”:
使用随机插值分解。识别矩阵转置的未在完全秩插值分解中使用的列。
None:
如果矩阵接近满秩,即矩阵秩与行数之差小于五,则使用“svd”。否则使用“pivot”。此默认行为可能会在未经通知的情况下更改。
默认值:无。对于输入稀疏的问题,此选项将被忽略,并使用基于“pivot”的算法,该算法见[5]。
对于特定方法的选项,请参见show_options('linprog')
。
x01-D 数组,可选
决策变量的猜测值,将由优化算法细化。当前仅由‘revised simplex’方法使用,并且仅当x0表示基本可行解时才能使用。
integrality1-D 数组或整数,可选
指示每个决策变量整数约束类型。
0
:连续变量;没有整数约束。
1
:整数变量;决策变量必须在bounds内为整数。
2
:半连续变量;决策变量必须在bounds内或取值0
。
3
:半整数变量;决策变量必须在bounds内为整数或取值0
。
默认情况下,所有变量均为连续的。
对于混合整数约束,请提供一个形状为c.shape的数组。为了从较短的输入推断出每个决策变量的约束条件,参数将使用np.broadcast_to广播到c.shape。
当前仅由'highs'
方法使用,否则忽略。
返回:
resOptimizeResult
一个scipy.optimize.OptimizeResult
,包含以下字段。请注意,字段的返回类型可能取决于优化是否成功,因此建议在依赖其他字段之前检查OptimizeResult.status:
x1-D 数组
使目标函数最小化同时满足约束条件的决策变量值。
funfloat
目标函数的最优值c @ x
。
slack1-D 数组
松弛变量(通常为正值),b_ub - A_ub @ x
。
con1-D 数组
等式约束的(通常为零的)残差,b_eq - A_eq @ x
。
successbool
找到最优解时为True
。
状态整数
表示算法退出状态的整数。
0
:优化成功终止。
1
:达到迭代限制。
2
:问题似乎无解。
3
:问题似乎无界。
4
:遇到数值困难。
nitint
所有阶段中执行的总迭代次数。
messagestr
表示算法退出状态的字符串描述。
另请参阅
show_options
求解器接受的附加选项。
注意
本节描述可以通过‘method’参数选择的可用求解器。
‘highs-ds’ 和 ‘highs-ipm’ 是 HiGHS 单纯形和内点法求解器的接口[13],‘highs’(默认)会自动选择两者之一。这些是 SciPy 中最快的线性规划求解器,特别适用于大型稀疏问题;哪个更快取决于问题本身。其他求解器(‘interior-point’、‘revised simplex’ 和 ‘simplex’)是遗留方法,将在 SciPy 1.11.0 中移除。
highs-ds 方法是 C++ 高性能双修订单纯形实现(HSOL)的包装器[13],[14]。highs-ipm 方法是 C++ 内点法实现的包装器[13],它具有交叉路由,因此与单纯形求解器一样精确。highs 方法会自动选择两者中的一种。对于涉及 linprog
的新代码,建议明确选择这三种方法值之一。
自版本 1.6.0 新增。
interior-point 方法使用在[4]中概述的原始-对偶路径跟踪算法。此算法支持稀疏约束矩阵,对于大型稀疏问题特别快速。然而,返回的解可能比单纯形方法稍微不准确,并且通常不与约束定义的多面体顶点对应。
自版本 1.0.0 新增。
修订单纯形法 方法使用修订的单纯形法,如[9]中所述,但在算法的每次迭代中,使用基础矩阵的因子分解[11]来有效地维护和解决线性系统。
自版本 1.3.0 新增。
simplex 方法使用 Dantzig 单纯形算法的传统全表实现[1],[2](不是 Nelder-Mead 单纯形)。此算法包含以保持向后兼容性和教育目的。
自版本 0.15.0 新增。
在应用 interior-point、revised simplex 或 simplex 之前,基于[8]的预处理过程尝试识别平凡的不可行性、平凡的无界性和潜在的问题简化。具体来说,它检查以下情况:
-
A_eq
或A_ub
中的零行,表示平凡约束; -
A_eq
和A_ub
中的零列,表示无约束变量; -
列单体在
A_eq
中,表示固定变量; -
列单体在
A_ub
中,表示简单边界。
如果预处理显示问题无界(例如,一个无约束和无界变量具有负成本)或不可行(例如,A_eq
中的零行与b_eq
中的非零对应),求解器将以适当的状态代码终止。请注意,预处理一旦检测到任何无界性的迹象就会终止;因此,当实际上问题是不可行时(但尚未检测到不可行性),可能会报告问题是无界的。因此,如果重要的是知道问题实际上是否不可行,请使用选项presolve=False
重新解决问题。
如果在预处理的单次通行中既未检测到不可行性也未检测到无界性,则在可能的情况下收紧界限并从问题中删除固定变量。然后,删除A_eq
矩阵的线性相关行(除非它们代表不可行性),以避免主要求解例程中的数值困难。请注意,几乎线性相关的行(在规定的容差内)也可以被删除,这在极少数情况下可能会改变最优解。如果这是一个问题,请从您的问题表达中消除冗余并使用选项rr=False
或presolve=False
运行。
这里可以进行几个潜在的改进:应该实现在[8]中概述的额外预处理检查,应该多次运行预处理例程(直到无法进一步简化为止),并且应该在冗余删除程序中实现[5]的更多效率改进。
经过预处理后,通过将(加紧的)简单界限转换为上界约束,为不等式约束引入非负松弛变量,并将无界变量表示为两个非负变量的差异,问题转换为标准形式。可选地,问题通过均衡转换自动进行缩放[12]。所选算法解决标准形式问题,并通过后处理程序将结果转换为原问题的解决方案。
参考文献
[1]
Dantzig, George B., 线性规划及其扩展。兰德公司研究学习普林斯顿大学出版社,普林斯顿,新泽西州,1963 年。
[2]
Hillier, S.H. 和 Lieberman, G.J. (1995), “数学规划导论”, 麦格劳-希尔, 第四章。
[3]
Bland, Robert G. 简单法的新有限枢轴规则。运筹学数学(2),1977 年:103-107 页。
[4]
Andersen, Erling D. 和 Knud D. Andersen. “MOSEK 内点优化器用于线性规划:同质算法的实现”。高性能优化。斯普林格美国出版,2000 年。197-232 页。
[5] (1,2,3)
Andersen, Erling D. “在大规模线性规划中找到所有线性相关的行。” 优化方法和软件 6.3(1995):219-227 页。
[6]
弗洛伊德,罗伯特·M。“基于牛顿方法的线性规划原始-对偶内点方法。”未发表的课程笔记,2004 年 3 月。可在 2017 年 2 月 25 日访问 ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf
[7]
福勒,罗伯特。“通过内点方法解线性规划问题。”未发表的课程笔记,2005 年 8 月 26 日。可在 2017 年 2 月 25 日访问 www.4er.org/CourseNotes/Book%20B/B-III.pdf
[8] (1,2)
安德森,埃尔林·D.,和克努德·D. 安德森。“线性规划中的预处理。”数学规划 71.2 (1995): 221-245。
[9]
贝茨马斯,迪米特里斯,和 J. Tsitsiklis。“线性规划导论。”Athena Scientific 1 (1997): 997。
[10]
安德森,埃尔林·D.,等人。大规模线性规划内点方法的实现。HEC/日内瓦大学,1996 年。
[11]
巴特尔斯,理查德·H。“单纯形方法的稳定化。”Journal in Numerische Mathematik 16.5 (1971): 414-434。
[12]
汤姆林,J. A.。“关于缩放线性规划问题。”数学规划研究 4 (1975): 146-166。
[13] (1,2,3)
黄甫,秦,加拉博娃,伊娃,费尔德迈尔,马克,和霍尔,J. A. J.。“HiGHS - 用于线性优化的高性能软件。”highs.dev/
[14]
黄甫,秦,和霍尔,J. A. J.。“对偶修订单纯形方法的并行化。”数学规划计算,10 (1),119-142,2018 年。DOI: 10.1007/s12532-017-0130-5
示例
考虑以下问题:
[\begin{split}\min_{x_0, x_1} \ -x_0 + 4x_1 & \ \mbox{这样的} \ -3x_0 + x_1 & \leq 6,\ -x_0 - 2x_1 & \geq -4,\ x_1 & \geq -3.\end{split}]
问题的表述形式不符合linprog
接受的格式。通过将“大于”不等式约束转换为“小于”不等式约束,同时将两边乘以(-1)的方法可以轻松解决这个问题。还需注意,最后一个约束实际上是简单的界限条件(-3 \leq x_1 \leq \infty)。最后,由于(x_0)没有边界,我们必须明确指定边界(-\infty \leq x_0 \leq \infty),因为默认情况下变量是非负的。将系数收集到数组和元组中后,该问题的输入为:
>>> from scipy.optimize import linprog
>>> c = [-1, 4]
>>> A = [[-3, 1], [1, 2]]
>>> b = [6, 4]
>>> x0_bounds = (None, None)
>>> x1_bounds = (-3, None)
>>> res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])
>>> res.fun
-22.0
>>> res.x
array([10., -3.])
>>> res.message
'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
边际(也称为对偶值 / 影子价格 / 拉格朗日乘子)和剩余量(余量)也是可用的。
>>> res.ineqlin
residual: [ 3.900e+01 0.000e+00]
marginals: [-0.000e+00 -1.000e+00]
例如,因为与第二个不等式约束相关联的边际为 -1,我们预期如果我们在第二个不等式约束的右侧添加一个小量eps
,目标函数的最优值将减少eps
:
>>> eps = 0.05
>>> b[1] += eps
>>> linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds]).fun
-22.05
Also, because the residual on the first inequality constraint is 39, we can decrease the right hand side of the first constraint by 39 without affecting the optimal solution.
>>> b = [6, 4] # reset to original values
>>> b[0] -= 39
>>> linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds]).fun
-22.0