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=None, jac_sparsity=None, max_nfev=None, verbose=0, args=(), kwargs=None)[源代码]#
求解带有变量约束的非线性最小二乘问题。
给定残差 f(x)（n 个实变量的 m 维实函数）和损失函数 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 维类数组或标量。如果参数x是复数或函数fun返回复数残差，则必须将其包装在实数参数的实数函数中，如“示例”部分末尾所示。- x0形状为 (n,) 或浮点数的类数组
自变量的初始猜测。如果为浮点数，则将其视为具有一个元素的一维数组。当 method 为 ‘trf’ 时，初始猜测可能会略微调整以使其充分位于给定的 bounds 内。
- jac{‘2-point’, ‘3-point’, ‘cs’, 可调用对象}, 可选
计算雅可比矩阵的方法（一个 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 元组或
Bounds，可选 有两种指定边界的方法
Bounds类的实例自变量的下限和上限。默认为无边界。每个数组必须与 x0 的大小匹配，或者为标量，在后一种情况下，所有变量的边界都相同。使用带有适当符号的
np.inf来禁用所有或某些变量的边界。
- method{‘trf’, ‘dogbox’, ‘lm’}, 可选
执行最小化的算法。
‘trf’ : 信赖域反射算法，特别适用于具有边界的大型稀疏问题。通常是稳健的方法。
‘dogbox’ : 带有矩形信赖域的 dogleg 算法，典型用例是带有边界的小型问题。不建议用于雅可比矩阵秩亏的问题。
‘lm’ : MINPACK 中实现的 Levenberg-Marquardt 算法。不处理边界和稀疏雅可比矩阵。通常是小型无约束问题最有效的方法。
默认为 ‘trf’。有关更多信息，请参阅“注释”。
- ftol浮点数或 None，可选
通过成本函数的变化来终止的容差。默认为 1e-8。当
dF < ftol * F时，且在最后一步中局部二次模型与真实模型之间有足够的一致性，则停止优化过程。如果为 None 且 ‘method’ 不是 ‘lm’，则禁用此条件终止。如果 ‘method’ 为 ‘lm’，则此容差必须高于机器精度。
- xtol浮点数或 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’，则此容差必须高于机器精度。
- gtol浮点数或 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’，则此容差必须高于机器精度。
- x_scale类数组或 ‘jac’，可选
每个变量的特征尺度。设置 x_scale 等同于在缩放变量中重新表述问题
xs = x / x_scale。另一种观点是，沿第 j 个维度的信赖域大小与x_scale[j]成正比。通过设置 x_scale，使得沿任何缩放变量的给定大小的步长对成本函数具有相似的影响，可以实现更好的收敛性。如果设置为 ‘jac’，则使用雅可比矩阵的列的逆范数迭代更新尺度（如 [JJMore] 中所述）。- lossstr 或 可调用对象，可选
确定损失函数。允许使用以下关键字值：
‘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’ 相似的属性。
如果为可调用对象，则它必须接受一个 1 维 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 将被视为用于有限差分方案的机器 epsilon 的常规“最佳”幂 [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’ 选项（bool，默认为 True），它向法方程添加正则化项，如果雅可比矩阵是秩亏的，则可以改善收敛性 [Byrd] (eq. 3.4)。
- jac_sparsity{None, 类数组对象, 稀疏矩阵}，可选
定义用于有限差分估计的雅可比矩阵的稀疏结构，其形状必须为 (m, n)。如果雅可比矩阵在每行中只有少量非零元素，则提供稀疏结构将大大加快计算速度 [Curtis]。零条目表示雅可比矩阵中的相应元素完全为零。如果提供，则强制使用 ‘lsmr’ 信赖域求解器。如果为 None（默认），则将使用密集差分。对 ‘lm’ 方法无效。
- verbose{0, 1, 2}，可选
算法的详细程度级别
0（默认）：静默工作。
1：显示终止报告。
2：显示迭代过程中的进度（‘lm’ 方法不支持）。
- args, kwargstuple 和 dict，可选
传递给 fun 和 jac 的其他参数。默认情况下都为空。调用签名是
fun(x, *args, **kwargs)，jac 也相同。
- 返回：
- resultOptimizeResult
定义了以下字段的
OptimizeResult- xndarray，形状 (n,)
找到的解。
- costfloat
解的成本函数值。
- funndarray，形状 (m,)
解处的残差向量。
- jacndarray，稀疏矩阵或 LinearOperator，形状 (m, n)
解处修改的雅可比矩阵，因为 J^T J 是成本函数黑塞矩阵的高斯-牛顿近似。类型与算法使用的类型相同。
- gradndarray，形状 (m,)
解处的成本函数梯度。
- optimalityfloat
一阶最优性度量。在无约束问题中，它始终是梯度的均匀范数。在约束问题中，它是与迭代期间的 gtol 进行比较的量。
- active_maskint 的 ndarray，形状 (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
如果满足其中一个收敛标准（status > 0），则为 True。
注释
方法 ‘lm’ (Levenberg-Marquardt) 调用 MINPACK (lmder, lmdif) 中实现的最小二乘算法的包装器。它运行公式化为信赖域类型算法的 Levenberg-Marquardt 算法。该实现基于论文 [JJMore]，它非常稳健，并且具有许多巧妙的技巧。对于无约束问题，它应该是您的首选。请注意，它不支持边界。此外，当 m < n 时它不起作用。
方法 ‘trf’（信赖域反射）的动机是解决方程组的过程，该方程组构成了如 [STIR] 中所述的带约束的最小化问题的一阶最优性条件。该算法迭代地解决信赖域子问题，并通过特殊的对角二次项进行增强，并且信赖域的形状由与边界的距离和梯度的方向决定。这种增强有助于避免直接步入边界，并有效地探索变量的整个空间。为了进一步提高收敛性，该算法考虑从边界反射的搜索方向。为了遵守理论要求，该算法保持迭代严格可行。对于稠密雅可比矩阵，信赖域子问题通过一种与 [JJMore] 中描述的方法非常相似的精确方法解决（并在 MINPACK 中实现）。与 MINPACK 实现的不同之处在于，雅可比矩阵的奇异值分解在每次迭代中执行一次，而不是 QR 分解和一系列 Givens 旋转消除。对于大型稀疏雅可比矩阵，使用了求解信赖域子问题的二维子空间方法 [STIR], [Byrd]。该子空间由缩放的梯度和由
scipy.sparse.linalg.lsmr提供的近似高斯-牛顿解构成。当没有施加约束时，该算法与 MINPACK 非常相似，并且通常具有可比的性能。该算法在无界和有界问题中都非常稳健，因此被选为默认算法。方法 ‘dogbox’ 在信赖域框架中运行，但考虑矩形信赖域，而不是传统的椭圆体 [Voglis]。当前信赖域和初始边界的交集仍然是矩形的，因此在每次迭代中，使用 Powell 的 dogleg 方法 [NumOpt] 近似求解受边界约束的二次最小化问题。对于稠密雅可比矩阵，可以精确计算所需的高斯-牛顿步长，或者对于大型稀疏雅可比矩阵，可以通过
scipy.sparse.linalg.lsmr近似计算。当雅可比矩阵的秩小于变量数时，该算法可能表现出缓慢的收敛速度。在变量数量较少的有界问题中，该算法通常优于 ‘trf’。鲁棒损失函数按照 [BA] 中的描述实现。其思想是在每次迭代时修改残差向量和雅可比矩阵，使得计算出的梯度和高斯-牛顿黑塞矩阵近似匹配成本函数的真实梯度和黑塞矩阵近似。然后该算法以正常方式进行，即，鲁棒损失函数被实现为标准最小二乘算法的简单包装。
在版本 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 et. al., “Numerical Recipes. The Art of Scientific Computing. 3rd edition”, Sec. 5.7.
[Byrd] (1,2)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 参数以bounds=([-np.inf, 1.5], np.inf)的形式指定给least_squares。我们还提供了分析雅可比矩阵
>>> 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])
因此，我们优化了 2n 个实变量的 2m 维实函数，而不是原始的 n 个复变量的 m 维复函数
>>> 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)