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=None, 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, callback=None, workers=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) 的目的是减少离群值（outliers）对解的影响。
- 参数:
- fun可调用对象
计算残差向量的函数，签名格式为
fun(x, *args, **kwargs)，即最小化过程针对其第一个参数进行。传递给该函数的参数x是一个形状为 (n,) 的 ndarray（即使 n=1 也绝不是标量）。它必须分配并返回一个形状为 (m,) 的类数组（array_like）或一个标量。如果参数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 能正确处理复数输入且可以在复平面上解析延拓的情况。如果是可调用对象，则按
jac(x, *args, **kwargs)调用，并应返回雅可比矩阵的良好近似（或精确值），其形式可以是类数组（应用 np.atleast_2d）、稀疏数组（为了性能推荐 csr_array）或scipy.sparse.linalg.LinearOperator。1.16.0 版本中的更改：增加了在 ‘lm’ 方法中使用 ‘3-point’、‘cs’ 关键字的能力。以前 ‘lm’ 仅限于 ‘2-point’ 和可调用对象。
- bounds类数组的 2 元组或
Bounds, 可选 指定边界有两种方式：
Bounds类的实例独立变量的上下限。默认无限制。每个数组必须与 x0 的大小匹配，或者是标量（此时所有变量的边界相同）。使用带有适当符号的
np.inf可以禁用所有或部分变量的边界。
- method{‘trf’, ‘dogbox’, ‘lm’}, 可选
执行最小化的算法。
‘trf’：信赖域反射算法（Trust Region Reflective），特别适用于带边界的大型稀疏问题。通常是一种稳健的方法。
‘dogbox’：具有矩形信赖域的 dogleg 算法，典型用例是带边界的小型问题。不推荐用于雅可比矩阵秩亏的问题。
‘lm’：MINPACK 中实现的 Levenberg-Marquardt 算法。不支持边界和稀疏雅可比矩阵。通常是处理小型无约束问题最有效的方法。
默认为 ‘trf’。更多信息请参见“注意”部分。
- ftol浮点数或 None, 可选
通过成本函数的变化来终止的容差。默认值为 1e-8。当满足
dF < ftol * F，且最后一步中局部二次模型与真实模型之间有足够的契合度时，优化过程停止。如果为 None 且 ‘method’ 不是 ‘lm’，则禁用此条件的终止。如果 ‘method’ 为 ‘lm’，则此容差必须高于机器 epsilon。
- 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’，则此容差必须高于机器 epsilon。
- 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’，则此容差必须高于机器 epsilon。
- x_scale{None, array_like, ‘jac’}, 可选
每个变量的特征尺度。设置 x_scale 相当于将问题重新表述为缩放变量
xs = x / x_scale。另一种观点是，信赖域沿第 j 维的大小与x_scale[j]成比例。通过设置 x_scale 使得沿任何缩放变量的给定大小的步骤对成本函数产生相似的影响，可以提高收敛性。如果设置为 ‘jac’，则使用雅可比矩阵列范数的倒数迭代更新尺度（如 [JJMore] 中所述）。每种方法的默认缩放（即如果x_scale 为 None）如下对于 ‘trf’：
x_scale == 1对于 ‘dogbox’：
x_scale == 1对于 ‘lm’：
x_scale == 'jac'
1.16.0 版本中的更改：默认关键字值从 1 更改为 None，表示使用默认缩放方法。对于 ‘lm’ 方法，默认缩放从 1 更改为 ‘jac’。发现这能提供更好的性能，并且与
leastsq执行的缩放相同。- loss字符串或可调用对象, 可选
确定损失函数。允许使用以下关键字值：
‘linear’（默认）：
rho(z) = z。给出标准的最小二乘问题。‘soft_l1’：
rho(z) = 2 * ((1 + z)**0.5 - 1)。L1（绝对值）损失的光滑近似。通常是稳健最小二乘的良好选择。‘huber’：
rho(z) = z 如果 z <= 1 否则 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_scale浮点数, 可选
内点和离群残差之间的软边界值，默认值为 1.0。损失函数评估如下：
rho_(f**2) = C**2 * rho(f**2 / C**2)，其中C是 f_scale，而rho由 loss 参数决定。当loss='linear'时此参数无效，但对于其他 loss 值，它至关重要。- max_nfevNone 或整数, 可选
对于所有方法，此参数控制每种方法使用的最大函数评估次数，不包括数值近似雅可比矩阵所使用的次数。如果为 None（默认值），则自动选择值为 100 * n。
1.16.0 版本中的更改：对于 ‘lm’ 方法，默认值更改为 100 * n，无论雅可比矩阵是可调用的还是数值估计的。以前使用估计雅可比矩阵时的默认值是 100 * n * (n + 1)，因为该方法包含了估计中使用的评估次数。
- 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_options字典, 可选
传递给信赖域解算器的关键字选项。
tr_solver='exact'：忽略 tr_options。tr_solver='lsmr'：scipy.sparse.linalg.lsmr的选项。此外，method='trf'支持 ‘regularize’ 选项（布尔值，默认为 True），它向正规方程添加正则化项，如果雅可比矩阵秩亏，则可提高收敛性 [Byrd] (eq. 3.4)。
- jac_sparsity{None, 类数组, 稀疏数组}, 可选
定义雅可比矩阵用于有限差分估计的稀疏结构，其形状必须为 (m, n)。如果雅可比矩阵在*每一*行中只有极少数非零元素，提供稀疏结构将极大地加快计算速度 [Curtis]。零项意味着雅可比矩阵中的相应元素恒为零。如果提供，将强制使用 ‘lsmr’ 信赖域解算器。如果为 None（默认值），则将使用稠密差分。对 ‘lm’ 方法没有影响。
- verbose{0, 1, 2}, 可选
算法的冗余级别
0（默认）：静默运行。
1：显示终止报告。
2：显示迭代过程中的进度（‘lm’ 方法不支持）。
- args, kwargs元组和字典, 可选
传递给 fun 和 jac 的附加参数。默认均为空。调用签名为
fun(x, *args, **kwargs)，jac 相同。- callbackNone 或可调用对象, 可选
算法在每次迭代时调用的回调函数。这可用于在每一步打印或绘制优化结果，并根据某些用户定义的条件停止优化算法。仅为 trf 和 dogbox 方法实现。
签名格式为
callback(intermediate_result: OptimizeResult)intermediate_result 是一个 `scipy.optimize.OptimizeResult`，其中包含当前迭代优化的中间结果。
回调还支持如下签名：
callback(x)使用内省来确定调用哪种签名。
如果 callback 函数引发 StopIteration，优化算法将停止并返回状态码 -2。
1.16.0 版本中添加。
- workers类 map 的可调用对象, 可选
类 map 的可调用对象，例如 multiprocessing.Pool.map，用于并行评估数值微分。该评估以
workers(fun, iterable)的形式执行。1.16.0 版本中添加。
- 返回:
- resultOptimizeResult
OptimizeResult具有定义的以下字段：- xndarray, 形状 (n,)
找到的解。
- cost浮点数
解处的成本函数值。
- funndarray, 形状 (m,)
解处的残差向量。
- jacndarray、稀疏数组或 LinearOperator, 形状 (m, n)
解处修改后的雅可比矩阵，在此意义上 J^T J 是成本函数 Hessian 矩阵的 Gauss-Newton 近似。类型与算法所用的类型相同。
- gradndarray, 形状 (m,)
解处成本函数的梯度。
- optimality浮点数
一阶最优性度量。在无约束问题中，它始终是梯度的无穷范数。在约束问题中，它是迭代过程中与 gtol 进行比较的量。
- active_mask整数型 ndarray, 形状 (n,)
每个分量显示相应约束是否处于激活状态（即变量是否在边界上）：
0：约束未激活。
-1：下限激活。
1：上限激活。
对于 ‘trf’ 方法可能有些武断，因为它生成一系列严格可行的迭代，而 active_mask 是在容差阈值内确定的。
- nfev整数
完成的函数评估次数。此数字不包括用于数值雅可比近似的函数调用。
1.16.0 版本中的更改：对于 ‘lm’ 方法，不再包括数值雅可比近似中使用的函数调用次数。这是为了使所有方法保持一致。
- njev整数或 None
完成的雅可比评估次数。如果在 ‘lm’ 方法中使用数值雅可比近似，则将其设置为 None。
- status整数
算法终止的原因：
-2：由于回调引发 StopIteration 而终止。
-1：MINPACK 返回了不正确的输入参数状态。
0：超过了最大函数评估次数。
1：满足了 gtol 终止条件。
2：满足了 ftol 终止条件。
3：满足了 xtol 终止条件。
4：同时满足了 ftol 和 xtol 终止条件。
- message字符串
终止原因的口头描述。
- success布尔值
如果满足其中一个收敛标准（status > 0），则为 True。
附注
‘lm’（Levenberg-Marquardt）方法调用 MINPACK 中实现的最小二乘算法 (lmder) 的包装器。它运行表述为信赖域型算法的 Levenberg-Marquardt 算法。该实现基于论文 [JJMore]，非常稳健且高效，带有许多巧妙的技巧。对于无约束问题，它应该是你的首选。请注意它不支持边界。此外，当 m < n 时它也不起作用。
‘trf’（信赖域反射）方法的动机是求解一个方程组的过程，该方程组构成了 [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] 中所述。其思想是在每次迭代中修改残差向量和雅可比矩阵，使得计算出的梯度和 Gauss-Newton Hessian 近似与成本函数的真实梯度和 Hessian 近似相匹配。然后算法以正常方式进行，即稳健损失函数作为标准最小二乘算法的一个简单包装器来实现。
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]保持无约束。为此，我们向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_array >>> def sparsity_broyden(n): ... sparsity = lil_array((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 维复函数，而是优化 2n 个实变量的 2m 维实函数。
>>> 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)