SciPy 1.12 中文文档(二十四)
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{None, ‘exact’, ‘lsmr’},可选
在迭代过程中解决无界最小二乘问题的方法:
- ‘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},可选
算法详细程度:
- 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”, “iteration”},可选
找到函数的不动点的方法,默认为“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 教程。