当前位置: 首页>>代码示例>>用法及示例精选>>正文


python 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={})

解决带有变量边界的非线性least-squares问题。

给定残差f(x)(n个实变量的m-dimensional实函数)和损失函数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)的目的是减少异常值对解的影响。

参数:

funcallable

具有签名的计算残差矢量的函数fun(x, *args, **kwargs),即最小化会针对第一个参数进行。争论x传递给此函数的是形状为(n,)的ndarray(即使n = 1,也不要标量)。它必须返回形状(m,)或标量的1-d 数组。如果论点x是复杂的还是函数fun返回复杂的残差,必须将其包装在实参的实函数中,如示例部分结尾所示。

x0array_like with shape (n,) 或 float

对自变量的初步猜测。如果为float,则将其视为具有一个元素的一维数组。

jac{‘2-point’, ‘3-point’, ‘cs’, callable}, 可选参数

计算雅可比矩阵的方法(m-by-n矩阵,其中元素(i,j)是f [i]相对于x [j]的偏导数)。关键字选择用于数字估计的有限差分方案。方案“ 3点”更准确,但所需的操作量是“ 2点”(默认值)的两倍。方案‘cs’使用复杂的步骤,虽然可能是最精确的步骤,但仅在fun正确地处理了复杂的输入并且可以解析地继续到复杂的平面时才适用。方法‘lm’始终使用“两点”方案。如果可调用,则用作jac(x, *args, **kwargs)并应以数组(应用np.atleast_2d),稀疏矩阵或a的形式返回雅可比行列式的良好近似值(或精确值)。scipy.sparse.linalg.LinearOperator

bounds2-tuple of array_like, 可选参数

自变量的上下限。默认为无界。每个数组必须匹配x0的大小或为标量,在后一种情况下,所有变量的界限都相同。使用np.inf带有适当的符号以禁用所有或某些变量的界限。

method{‘trf’, ‘dogbox’, ‘lm’}, 可选参数

执行最小化的算法。

  • ‘trf’:Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method.

  • ‘dogbox’:dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian.

  • ‘lm’:Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t handle bounds and sparse Jacobians. Usually the most efficient method for small unconstrained problems.

默认值为‘trf’。有关更多信息,请参见注释。

ftolfloat 或 None, 可选参数

成本函数变化导致的终止公差。默认值为1e-8。优化过程在以下时间停止dF < ftol * F,并且在最后一步中,本地二次模型与真实模型之间有足够的一致性。如果为None,则禁用此条件的终止。

xtolfloat 或 None, 可选参数

通过改变自变量终止的公差。默认值为1e-8。确切条件取决于所使用的方法:

  • For ‘trf’ and ‘dogbox’:norm(dx) < xtol * (xtol + norm(x))

  • For ‘lm’:Delta < xtol * norm(xs), where Delta is a trust-region radius and xs is the value of x scaled according to x_scale parameter (see below).

如果为None,则禁用此条件的终止。

gtolfloat 或 None, 可选参数

渐变规范终止的公差。默认值为1e-8。确切条件取决于所使用的方法:

  • For ‘trf’:norm(g_scaled, ord=np.inf) < gtol, where g_scaled is the value of the gradient scaled to account for the presence of the bounds [STIR].

  • For ‘dogbox’:norm(g_free, ord=np.inf) < gtol, where g_free is the gradient with respect to the variables which are not in the optimal state on the boundary.

  • For ‘lm’:the maximum absolute value of the cosine of angles between columns of the Jacobian and the residual vector is less than gtol, or the residual vector is zero.

如果为None,则禁用此条件的终止。

x_scalearray_like 或 ‘jac’, 可选参数

每个变量的特征量表。设置x_scale等效于在比例变量中重新构造问题xs = x / x_scale。另一种观点是,信任区域沿j-th维度的大小与x_scale[j]。通过设置可以提高收敛性x_scale这样,沿着任何比例变量的给定大小的步长对成本函数的影响就类似。如果设置为‘jac’,则使用雅可比矩阵列的反范数迭代缩放比例(如[JJ更多])。

lossstr 或 callable, 可选参数

确定损失函数。允许使用以下关键字值:

  • ‘linear’ (default):rho(z) = z. Gives a standard least-squares problem.

  • ‘soft_l1’:rho(z) = 2 * ((1 + z)**0.5 - 1). The smooth approximation of l1 (absolute value) loss. Usually a good choice for robust least squares.

  • ‘huber’:rho(z) = z if z <= 1 else 2*z**0.5 - 1. Works similarly to ‘soft_l1’.

  • ‘cauchy’:rho(z) = ln(1 + z). Severely weakens outliers influence, but may cause difficulties in optimization process.

  • ‘arctan’:rho(z) = arctan(z). Limits a maximum loss on a single residual, has properties similar to ‘cauchy’.

如果可调用,则必须采用一维ndarrayz=f**2并返回形状(3,m)的数组,其中第0行包含函数值,第1行包含一阶导数,第2行包含二阶导数。方法‘lm’仅支持‘linear’丢失。

f_scalefloat, 可选参数

异常值和异常值之间的软裕度值,默认值为1.0。损失函数的评估如下rho_(f**2) = C**2 * rho(f**2 / C**2),在其中C是f_scale,并且rho由...决定失利参数。此参数对loss='linear',但其他失利重视它至关重要。

max_nfevNone 或 int, 可选参数

终止之前的最大函数评估数。如果为None(默认),则自动选择该值:

  • For ‘trf’ and ‘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。如果为无(默认值),则对于使用的有限差分方案,diff_step被视为机器epsilon的常规“optimal”幂[NR]

tr_solver{None, ‘exact’, ‘lsmr’}, 可选参数

解决trust-region个子问题的方法,仅与‘trf’和‘dogbox’方法有关。

  • ‘exact’ is suitable for not very large problems with dense Jacobian matrices. The computational complexity per iteration is comparable to a singular value decomposition of the Jacobian matrix.

  • ‘lsmr’ is suitable for problems with sparse and large Jacobian matrices. It uses the iterative procedure scipy.sparse.linalg.lsmr for finding a solution of a linear least-squares problem and only requires matrix-vector product evaluations.

如果为None(默认),则根据第一次迭代返回的Jacobian类型选择求解器。

tr_optionsdict, 可选参数

关键字选项已传递给trust-region解算器。

  • tr_solver='exact':tr_options are ignored.

  • tr_solver='lsmr':options for scipy.sparse.linalg.lsmr. Additionally method='trf' supports ‘regularize’ option (bool, default is True) which adds a regularization term to the normal equation, which improves convergence if the Jacobian is rank-deficient [Byrd] (eq. 3.4).

jac_sparsity{None, array_like, sparse matrix}, 可选参数

定义用于有限差分估计的Jacobian矩阵的稀疏结构,其形状必须为(m,n)。如果雅可比行列每行只有几个非零元素,则提供稀疏结构将大大加快计算速度[柯蒂斯]。零条目表示雅可比行列中的对应元素相同为零。如果提供,则强制使用‘lsmr’ trust-region解算器。如果没有(默认),则将使用密集差分。对‘lm’方法无效。

verbose{0, 1, 2}, 可选参数

算法的详细程度:

  • 0 (default):work silently.

  • 1:display a termination report.

  • 2:display progress during iterations (not supported by ‘lm’ method).

args, kwargstuple and dict, 可选参数

其他参数传递给fun和jac。两者默认为空。主叫签名是fun(x, *args, **kwargs)和相同的江淮

返回值:

OptimizeResult定义了以下字段:
xndarray,形状(n,)

找到解决方案。

cost浮动

成本函数在解决方案中的价值。

funndarray,形状(m,)

溶液中的残差矢量。

jacndarray,稀疏矩阵或LinearOperator,形状(m,n)

在解决方案上,在J^T J是成本函数的黑森州的Gauss-Newton近似的意义上,修改了雅可比矩阵。类型与算法使用的类型相同。

gradndarray,形状(m,)

解决方案中成本函数的梯度。

optimality浮动

First-order最佳度量。在无约束的问题中,它始终是梯度的统一范数。在受约束的问题中,在迭代过程中将其与gtol进行比较。

active_maskint的ndarray,形状(n,)

每个组件都显示相应的约束是否处于活动状态(即,变量是否处于边界):

  • 0:a constraint is not active.

  • -1:a lower bound is active.

  • 1:an upper bound is active.

对于‘trf’方法,它可能在某种程度上是任意的,因为它生成严格可行的迭代序列,并且在公差阈值内确定active_mask。

nfev整型

完成函数评估的次数。与‘lm’方法相反,方法‘trf’和‘dogbox’不计算数字雅可比近似的函数调用。

njev整型或无

完成的雅可比评估数。如果在‘lm’方法中使用了数字Jacobian近似,则将其设置为None。

status整型

算法终止的原因:

  • -1:improper input parameters status returned from MINPACK.

  • 0:the maximum number of function evaluations is exceeded.

  • 1:gtol termination condition is satisfied.

  • 2:ftol termination condition is satisfied.

  • 3:xtol termination condition is satisfied.

  • 4:Both ftol and xtol termination conditions are satisfied.

message力量

终止原因的口头描述。

success布尔

如果满足收敛标准之一(状态> 0),则为True。

注意:

方法‘lm’(Levenberg-Marquardt)通过MINPACK(lmder,lmdif)中实现的least-squares算法调用包装器。它运行公式化为trust-region类型算法的Levenberg-Marquardt算法。实施是基于纸的[JJ更多],它具有许多巧妙的技巧,非常强大且高效。对于不受限制的问题,它应该是您的首选。请注意,它不支持边界。当m <n时,它也不起作用。

方法‘trf’(可信区域反射)是由求解方程组的过程引起的,这些方程组构成了bound-constrained最小化问题的first-order最优条件,公式如下:[搅拌]。该算法迭代地解决了trust-region子问题,该子问题由特殊的对角二次项扩展,并且具有trust-region形状,该形状由距边界的距离和梯度方向确定。此增强函数有助于避免直接步入界限,并有效地探索变量的整个空间。为了进一步提高收敛性,该算法考虑了从边界反映的搜索方向。为了遵守理论要求,该算法保持迭代严格可行。对于密集的Jacobian问题,trust-region子问题可以通过一种精确的方法来解决,该方法非常类似于[JJ更多](并在MINPACK中实现)。与MINPACK实现不同的是,雅各比矩阵的奇异值分解每次迭代执行一次,而不是QR分解和一 Series Givens旋转消除。对于较大的稀疏雅可比矩阵,使用求解trust-region子问题的二维子空间方法[搅拌][伯德]。子空间由缩放梯度和由以下对象提供的近似Gauss-Newton解决方案覆盖scipy.sparse.linalg.lsmr。当没有施加任何约束时,该算法与MINPACK非常相似,并且通常具有可比的性能。该算法在无界和有界问题上都非常健壮,因此被选为默认算法。

方法‘dogbox’在trust-region框架中运行,但是考虑了矩形信任区域,而不是常规的椭圆体[Voglis]。当前信任区域和初始边界的交集再次为矩形,因此在每次迭代中,都会通过Powell的dogleg方法近似求解受边界约束的二次最小化问题[NumOpt]。所需的Gauss-Newton步骤可以针对密集的Jacobian精确计算,或者近似计算为scipy.sparse.linalg.lsmr用于大型稀疏的雅各布主义者。当雅可比行列的秩小于变量数时,该算法可能会表现出较慢的收敛性。在有少量变量的有界问题中,该算法通常优于‘trf’。

稳健的损失函数如以下所述实现[BA]。想法是在每次迭代中修改残差矢量和Jacobian矩阵,以使计算出的梯度和Gauss-Newton Hessian近似值与成本函数的真实梯度和Hessian近似值匹配。然后,该算法以正常方式进行,即,将鲁棒损失函数实现为对标准least-squares算法的简单包装。

0.17.0版中的新函数。

参考文献:

STIR(123)

M. A. Branch,T。F. Coleman和Y. Li,“针对Large-Scale Bound-Constrained最小化问题的子空间,内部和共轭梯度方法,” SIAM科学计算杂志,第1卷。 21,第1号,第1-23页,1999年。

NR

威廉·H·普雷斯等等,“数字食谱。科学计算的艺术。第三版”,第二节5.7。

伯德(12)

R. H. Byrd,R。B. Schnabel和G. A. Shultz,“通过最小化二维子空间来近似解决信任区问题”,数学。 《程序设计》,第40卷,第247-263页,1988年。

柯蒂斯

A. Curtis,M。J. D. Powell和J. Reid,“关于稀疏Jacobian矩阵的估计”,数学与应用研究所学报,第13页,第117-120页,1974年。

JJ更多(123)

J. J. More,“ Levenberg-Marquardt算法:实现与理论”,《数值分析》,编辑。 G.A.Watson,《数学讲义630》(Springer Verlag),第105-116页,1977年。

沃格利斯

C. Voglis和I. E. Lagaris,“用于无约束和有界约束非线性优化的矩形信任区域Dogleg方法”,WSEAS应用数学国际会议,希腊科孚,2004年。

NumOpt

J. Nocedal和S. J. Wright,“数值优化,第二版”,第4章。

BA

B.特里格斯等等人,“捆绑调整-现代综合”,视觉算法国际研讨会论文集:理论与实践,第298-372页,1999年。

例子:

在此示例中,我们找到了Rosenbrock函数的最小值,该变量不受自变量的限制。

>>> 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.5x[0]不受约束。为此,我们将bounds参数指定为least_squares形式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三对角vector-valued函数求解方程组(即成本函数至少应为零):

>>> 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是要估计的参数。

首先,定义函数以生成带有噪声和异常值的数据,定义模型参数并生成数据:

>>> def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):
...     y = a + b * np.exp(t * c)
...
...     rnd = np.random.RandomState(random_state)
...     error = noise * rnd.randn(t.size)
...     outliers = rnd.randint(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])

计算标准的least-squares解决方案:

>>> 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))

最后绘制所有曲线。我们看到,通过选择适当的损失,即使存在强烈的异常值,我们也可以获得接近最佳值的估计值。但是请记住,通常建议首先尝试‘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()
../_images/scipy-optimize-least_squares-1_00_00.png

在下一个示例中,我们展示如何使用以下方法优化复杂变量的复值残差函数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-dimensional复杂函数,我们优化了2n个实变量的2m-dimensional实函数:

>>> 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.least_squares的API实现见:[源代码]

相关用法

注:本文由纯净天空筛选整理自 scipy.optimize.least_squares。非经特殊声明,原始代码版权归原作者所有,本译文的传播和使用请遵循“署名-相同方式共享 4.0 国际 (CC BY-SA 4.0)”协议。