scipy.optimize.

最小二乘法#

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)（一个 n 个实变量的 m-D 实函数）和损失函数 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(x, *args, **kwargs)，即最小化过程是相对于其第一个参数进行的。传递给此函数的参数 x 是一个形状为 (n,) 的 ndarray（即使 n=1 也不会是标量）。它必须分配并返回一个形状为 (m,) 的类数组对象或一个标量。如果参数 x 是复数或函数 fun 返回复数残差，则必须将其包装在一个实数函数的实数参数中，如示例部分末尾所示。

x0类似数组，形状为 (n,) 或浮点数

对独立变量的初始猜测。如果是浮点数，它将被视为一个包含一个元素的一维数组。当 method 为 ‘trf’ 时，初始猜测可能会稍微调整，以充分位于给定的 bounds 内。

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

计算雅可比矩阵（一个 m-by-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 : 数组类或 Bounds 的 2-tuple，可选数组类对象的2元组或

有两种指定边界的方法：

  1. Bounds 类的实例

  2. 独立变量的上下界。默认为无界。每个数组的大小必须与 x0 匹配，或者为标量，在后一种情况下，所有变量的界限将相同。使用 np.inf 并带有适当的符号来禁用所有或某些变量的界限。

方法{‘trf’, ‘dogbox’, ‘lm’}, 可选

执行最小化的算法。

  • ‘trf’ : 信赖域反射算法，特别适合带有边界的大型稀疏问题。通常是一种稳健的方法。

  • ‘dogbox’ : 带有矩形信任区域的狗腿算法，典型用例是带有边界的小问题。不推荐用于雅可比矩阵秩不足的问题。

  • ‘lm’ : 在 MINPACK 中实现的 Levenberg-Marquardt 算法。不处理边界和稀疏雅可比矩阵。通常是小型无约束问题的最有效方法。

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

ftol浮点数或无，可选

成本函数变化的终止容差。默认值为 1e-8。当 dF < ftol * F 时，优化过程停止，并且在最后一步中，局部二次模型与真实模型之间有足够的吻合。

如果为 None 且 ‘method’ 不是 ‘lm’，则此条件下的终止将被禁用。如果 ‘method’ 是 ‘lm’，则此容差必须大于机器 epsilon。

xtol浮点数或无，可选

对独立变量变化终止的容差。默认值为 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。

gtol浮点数或无，可选

梯度范数终止容差。默认值为 1e-8。确切的终止条件取决于所使用的 方法：

  • 对于 ‘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_scale类数组 或 ‘jac’, 可选

每个变量的特征尺度。设置 x_scale 等同于在缩放变量 xs = x / x_scale 中重新表述问题。另一种观点是，沿第 j 维的信任区域的大小与 x_scale[j] 成正比。通过设置 x_scale 使得沿任何缩放变量的大小步长对成本函数的影响相似，可以实现更好的收敛。如果设置为 ‘jac’，则使用雅可比矩阵列的逆范数迭代更新尺度（如 [JJMore] 中所述）。

损失str 或 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) 的类数组对象，其中第 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。

  • 对于 ‘lm’ : 如果 jac 是可调用的，则为 100 * n，否则为 100 * n * (n + 1)（因为 ‘lm’ 在雅可比估计中计算函数调用）。

diff_stepNone 或类数组，可选

确定用于雅可比矩阵有限差分近似的相对步长。实际步长计算为 x * diff_step。如果为 None（默认），则 diff_step 被视为所用有限差分方案的机器精度的常规“最优”幂次 [NR]。

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

解决信赖域子问题的方法，仅与 ‘trf’ 和 ‘dogbox’ 方法相关。

  • ‘exact’ 适用于具有稠密雅可比矩阵的非非常大问题。每次迭代的计算复杂度与雅可比矩阵的奇异值分解相当。

  • ‘lsmr’ 适用于具有稀疏和大雅可比矩阵的问题。它使用迭代过程 scipy.sparse.linalg.lsmr 来寻找线性最小二乘问题的解，并且只需要矩阵-向量乘积的评估。

如果为 None（默认），求解器将根据第一次迭代返回的雅可比矩阵类型来选择。

tr_optionsdict, 可选

传递给信赖域求解器的关键词选项。

  • tr_solver='exact': tr_options 被忽略。

  • tr_solver='lsmr': scipy.sparse.linalg.lsmr 的选项。此外，method='trf' 支持 ‘regularize’ 选项（布尔值，默认为 True），该选项向正规方程添加一个正则化项，如果雅可比矩阵是秩亏的，这可以改善收敛性 [R20fc1df64af7-Byrd]_（方程 3.4）。

jac_sparsity{None, 类数组, 稀疏矩阵}, 可选

定义有限差分估计的雅可比矩阵的稀疏结构，其形状必须为 (m, n)。如果雅可比矩阵在 每一行 中只有少数非零元素，提供稀疏结构将大大加快计算速度 [Curtis]。零条目意味着雅可比矩阵中相应元素为零。如果提供，将强制使用 ‘lsmr’ 信赖域求解器。如果为 None（默认），则将使用密集差分。对 ‘lm’ 方法无效。

详细{0, 1, 2}, 可选

算法详细程度的级别：

  • 0 (默认) : 静默工作。

  • 1 : 显示终止报告。

  • 2 : 在迭代过程中显示进度（’lm’ 方法不支持）。

args, kwargstuple 和 dict，可选

传递给 fun 和 jac 的额外参数。默认情况下均为空。调用签名是 fun(x, *args, **kwargs)，jac 也是如此。

返回:
结果优化结果

OptimizeResult 具有以下定义的字段：

xndarray, 形状 (n,)

找到解决方案。

成本浮动

解处的成本函数值。

乐趣ndarray, 形状 (m,)

解处的残差向量。

jacndarray、稀疏矩阵或线性算子，形状为 (m, n)

在解处的修正雅可比矩阵，即 J^T J 是成本函数 Hessian 的高斯-牛顿近似。类型与算法使用的类型相同。

梯度ndarray, 形状 (m,)

解处的成本函数的梯度。

最优性浮动

一阶最优性度量。在无约束问题中，它始终是梯度的均匀范数。在有约束问题中，它是迭代过程中与 gtol 进行比较的量。

active_mask整数 ndarray，形状为 (n,)

每个组件显示相应的约束是否处于活动状态（即变量是否处于边界）：

  • 0 : 约束不活跃。

  • -1 : 一个下界是活跃的。

  • 1 : 一个上界是活动的。

对于 ‘trf’ 方法来说，这可能有些随意，因为它生成了一系列严格可行的迭代，而 active_mask 是在一个容差阈值内确定的。

nfev整数

已完成的功能评估次数。’trf’ 和 ‘dogbox’ 方法不计算数值雅可比近似的函数调用次数，这与 ‘lm’ 方法相反。

njev整数或无

完成的雅可比矩阵评估次数。如果在 ‘lm’ 方法中使用了数值雅可比矩阵近似，则设置为 None。

状态整数

算法终止的原因：

  • -1 : MINPACK 返回了不正确的输入参数状态。

  • 0 : 函数评估的最大次数已超过。

  • 1 : 满足 gtol 终止条件。

  • 2 : ftol 终止条件已满足。

  • 3 : xtol 终止条件已满足。

  • 4 : ftol 和 xtol 终止条件都已满足。

消息str

终止原因的口头描述。

成功布尔

如果满足其中一个收敛标准则为真（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 旋转消除。对于大型稀疏雅可比矩阵，使用二维子空间方法求解信赖域子问题 [STIR]，[Byrd]。该子空间由缩放梯度和 scipy.sparse.linalg.lsmr 提供的近似 Gauss-Newton 解张成。当没有施加约束时，该算法与 MINPACK 非常相似，并且通常具有可比性能。该算法在无界和有界问题中都表现得相当稳健，因此被选为默认算法。

方法 ‘dogbox’ 在信赖域框架下操作，但考虑的是矩形信赖域，而不是传统的椭球形 [Voglis]。当前信赖域与初始边界的交集再次是矩形的，因此在每次迭代中，通过 Powell 的 dogleg 方法 [NumOpt] 近似求解受边界约束的二次最小化问题。所需的 Gauss-Newton 步长可以精确计算密集的雅可比矩阵，或者通过 scipy.sparse.linalg.lsmr 近似计算大型稀疏雅可比矩阵。当雅可比矩阵的秩小于变量数量时，该算法可能会表现出缓慢的收敛性。在变量数量较少的边界问题中，该算法通常优于 ‘trf’。

鲁棒损失函数按照 [BA] 中的描述实现。其思想是在每次迭代中修改残差向量和雅可比矩阵，使得计算出的梯度和高斯-牛顿海森近似与成本函数的真实梯度和海森近似相匹配。然后算法以正常方式进行，即，鲁棒损失函数作为标准最小二乘算法的简单包装来实现。

Added in version 0.17.0.

参考文献

[STIR] (1,2,3)

M. A. Branch, T. F. Coleman, and Y. Li, “A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems,” SIAM Journal on Scientific Computing, Vol. 21, Number 1, pp 1-23, 1999.

[NR]

William H. Press 等人，《数值秘笈：科学计算的艺术。第三版》，第5.7节。

[Byrd]

R. H. Byrd, R. B. Schnabel and 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, and 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, ed. G. A. Watson, Lecture Notes in Mathematics 630, Springer Verlag, pp. 105-116, 1977.

[Voglis]

C. Voglis and I. E. Lagaris, “A Rectangular Trust Region Dogleg Approach for Unconstrained and Bound Constrained Nonlinear Optimization”, WSEAS International Conference on Applied Mathematics, Corfu, Greece, 2004.

[NumOpt]

J. Nocedal and S. J. Wright, “Numerical optimization, 2nd edition”, Chapter 4.

[BA]

B. Triggs et. al., “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] 不受约束。为此，我们将 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三对角向量值函数：

>>> 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()
../../_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])

因此，我们优化的是一个 2m-D 实函数，而不是原始的 n 复变量的 m-D 复函数，这个实函数有 2n 个实变量：

>>> 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)
On this page