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={})[source]#

求解具有变量界限的非线性最小二乘问题。

给定残差 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,) 的类数组或浮点数

独立变量的初始猜测。如果为浮点数，则将被视为具有一个元素的 1 维数组。当 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, 可选

有两种方法可以指定边界

  1. Bounds 类的实例

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

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

执行最小化的算法。

  • ‘trf’ : 信赖域反射算法，特别适用于具有界限的大型稀疏问题。通常是稳健的方法。

  • ‘dogbox’ : 具有矩形信赖域的狗腿算法，典型用例是具有界限的小型问题。不建议用于雅可比矩阵秩亏的问题。

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

默认值为 '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] 中所述）。

loss字符串或可调用对象，可选

确定损失函数。以下关键字值是允许的

  • ‘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 或 array_like，可选

确定雅可比矩阵有限差分近似值的相对步长。实际步长计算为 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’ 选项（bool，默认值为 True），该选项向正规方程添加正则化项，如果雅可比矩阵秩亏损，则可以改善收敛性 [Byrd]（等式 3.4）。

jac_sparsity{None, array_like, sparse matrix}，可选

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

verbose{0, 1, 2}，可选

算法的详细程度

  • 0（默认）：静默运行。

  • 1：显示终止报告。

  • 2：显示迭代过程（‘lm’ 方法不支持）。

args, kwargs元组和字典，可选

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

返回：:
resultOptimizeResult

OptimizeResult，其中定义了以下字段

xndarray，形状 (n,)

找到的解。

costfloat

解处成本函数的值。

funndarray，形状 (m,)

解处的残差向量。

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

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

gradndarray，形状 (m,)

解处的成本函数的梯度。

optimalityfloat

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

active_maskint 类型的 ndarray，形状 (n,)

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

  • 0：约束不处于活动状态。

  • -1：下界处于活动状态。

  • 1：上界处于活动状态。

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

nfevint

执行的函数评估次数。与 ‘lm’ 方法不同，‘trf’ 和 ‘dogbox’ 方法不计算用于数值雅可比矩阵近似的函数调用次数。

njevint 或 None

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

statusint

算法终止的原因

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

  • 0：超过最大函数评估次数。

  • 1：满足 gtol 终止条件。

  • 2：满足 ftol 终止条件。

  • 3：满足 xtol 终止条件。

  • 4：同时满足 ftol 和 xtol 终止条件。

messagestr

终止原因的文字描述。

successbool

如果满足其中一个收敛标准（status > 0），则为 True。

另请参见

leastsq

MINPACK 实现的 Levenberg-Marquadt 算法的旧版包装器。

curve_fit

应用于曲线拟合问题的最小二乘法。

注释

方法 ‘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] 近似地求解受限于边界的二次最小化问题。对于稠密的雅可比矩阵，可以精确地计算出所需的 Gauss-Newton 步，对于大型稀疏的雅可比矩阵，可以通过 scipy.sparse.linalg.lsmr 近似地计算。当雅可比矩阵的秩小于变量个数时，该算法可能会表现出缓慢的收敛速度。在具有少量变量的边界问题中，该算法通常优于 ‘trf’ 方法。

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

在版本 0.17.0 中添加。

参考文献

[STIR] (1,2,3)

M. A. Branch, T. F. Coleman 和 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 等，“Numerical Recipes. The Art of Scientific Computing. 3rd edition”, Sec. 5.7。

[Byrd] (1,2)

R. H. Byrd, R. B. Schnabel 和 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 和 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 和 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 和 S. J. Wright，“Numerical optimization, 2nd edition”, Chapter 4。

[BA]

B. Triggs 等，“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])

因此，我们不再优化原始的 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)
本页内容