简介与构造函数
在scipy中，非线性最小二乘法的目的是找到一组函数，使得误差函数的平方和最小，可以表示为如下公式
arg min f i F ( x ) = 0.5 ∑ i = 0 m − 1 ρ ( f i ( x ) 2 ) , x ∈ [ L , R ] \argmin_{f_i} F(x) = 0.5\sum_{i=0}^{m-1}\rho(f_i(x)^2),\quad x\in[L,R] fiargminF(x)=0.5i=0∑m−1ρ(fi(x)2),x∈[L,R]
其中 ρ \rho ρ表示损失函数，可以理解为对 f i ( x ) f_i(x) fi(x)的一次预处理。
scipy.optimize中封装了非线性最小二乘法函数least_squares，其定义为
least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, f_scale, loss, jac_sparsity, max_nfev, verbose, args, kwargs)
其中，func和x0为必选参数，func为待求解函数，x0为函数输入的初值，这两者无默认值，为必须输入的参数。
bound为求解区间，默认
(
−
∞
,
∞
)
(-\infty,\infty)
(−∞,∞)，verbose为1时，会有终止输出，为2时会print更多的运算过程中的信息。此外下面几个参数用于控制误差，比较简单。
| 默认值 | 备注 | |
|---|---|---|
| ftol | 1 0 − 8 10^{-8} 10−8 | 函数容忍度 |
| xtol | 1 0 − 8 10^{-8} 10−8 | 自变量容忍度 |
| gtol | 1 0 − 8 10^{-8} 10−8 | 梯度容忍度 |
| x_scale | 1.0 | 变量的特征尺度 |
| f_scale | 1.0 | 残差边际值 |
loss为损失函数，就是上面公式中的
ρ
\rho
ρ，默认为linear，可选值包括
-
linear： ρ ( z ) = z \rho(z)=z ρ(z)=z，此为标准非线性最小二乘法 -
soft_l1： ρ ( z ) = 2 ( 1 + z − 1 ) \rho(z)=2(\sqrt{1 + z}-1) ρ(z)=2(1+z−1)， 相当于L1损失的平滑 -
huber： 当 z ⩽ 1 z\leqslant1 z⩽1时， ρ ( z ) = z \rho(z)=z ρ(z)=z，否则 ρ ( z ) = 2 z − 1 \rho(z)=2\sqrt{z}-1 ρ(z)=2z−1， 表现与soft_l1相似 -
cauchy： ρ ( z ) = ln ( 1 + z ) \rho(z) = \ln(1 + z) ρ(z)=ln(1+z) -
arctan： ρ ( z ) = arctan ( z ) \rho(z) = \arctan(z) ρ(z)=arctan(z)
迭代策略
上面的公式仅给出了算法的目的，但并未暴露其细节。关于如何找到最小值，则需要确定搜索最小值的方法，method为最小值搜索的方案，共有三种选项，默认为trf
-
trf：即TrustRegion Reflective，信赖域反射算法 -
dogbox：信赖域狗腿算法 -
lm：Levenberg-Marquardt算法
这三种方法都是信赖域方法的延申，信赖域的优化思想其实就是从单点的迭代变成了区间的迭代，由于本文的目的是介绍scipy中所封装好的非线性最小二乘函数，故而仅对其原理做简略的介绍。
对于优化问题 min x ∈ R n f ( x ) \min_{x\in R^n} f(x) minx∈Rnf(x)而言，定义当前点的邻域
Ω k = { x ∈ R n ∣ ∥ x − x k ∥ ⩽ r } \Omega_k=\{x\in R^n\big\vert \Vert x-x_k\Vert\leqslant r\} Ωk={x∈Rn ∥x−xk∥⩽r}
其中 r r r为置信半径，假设在这个邻域内，目标函数可以近似为线性或二次函数，则可通过二次模型得到区间中的极小值点 s k s_k sk。然后以这个极小值点为中心，继续优化信赖域所对应的区间。
记 s = x − x k s=x-x_k s=x−xk, 表示步长； g k = ∇ f ( x k ) , B k ≈ ∇ 2 f ( x k ) g_k=\nabla f(x_k), B_k\approx\nabla^2 f(x_k) gk=∇f(xk),Bk≈∇2f(xk)分别是函数的一阶导数和黑塞矩阵的近似，则对上述问题进行Taylor展开，可以得到一个二阶的近似模型
arg min q k ( s ) = f ( x k ) + g k T s + 1 2 s T B k s , ∥ s ∥ ⩽ r \argmin q^{k}(s)=f(x_k)+g_k^Ts+\frac{1}{2}s^TB_ks, \Vert s\Vert\leqslant r argminqk(s)=f(xk)+gkTs+21sTBks,∥s∥⩽r
B k B_k Bk被定义为黑塞矩阵的近似，当其表示为雅可比矩阵 J k T J k J_k^TJ_k JkTJk时，所得到的迭代方案便是著名的高斯牛顿法，而LM算法在在高斯牛顿的基础上，添加了一个阻尼因子，可记作 J k T J k + μ I J_k^TJ_k+\mu I JkTJk+μI。
狗腿法的鲍威尔提出的一种迭代方案，特点是新定义了下降比，即随着 s s s的不断前进，目标函数和模型函数都会发生变化，则可定义二者的变化率 r k = f ( x k ) − f ( x k + s ) q k ( 0 ) − q k ( s ) r_k=\frac{f(x_k)-f(x_k+s)}{q^k(0)-q^k(s)} rk=qk(0)−qk(s)f(xk)−f(xk+s)，根据这个值的变化，来调整信赖域半径 r r r的值。
以上就是信赖域方法的基本原理。
雅可比矩阵
在了解了信赖域方法之后，就会明白雅可比矩阵在数值求解时的重要作用，而如何计算雅可比矩阵，则是接下来需要考虑的问题。jac参数为计算雅可比矩阵的方法，主要提供了三种方案，分别是基于两点的2-point；基于三点的3-point；以及基于复数步长的cs。一般来说，三点的精度高于两点，但速度也慢一倍。
此外，可以输入自定义函数来计算雅可比矩阵。文章来源：https://www.toymoban.com/news/detail-699787.html
测试
最后，测试一下非线性最小二乘法文章来源地址https://www.toymoban.com/news/detail-699787.html
import numpy as np
from scipy.optimize import least_squares
def test(xs):
_sum = 0.0
for i in range(len(xs)):
_sum = _sum + (1-np.cos((xs[i]*i)/5)*(i+1))
return _sum
x0 = np.random.rand(5)
ret = least_squares(test, x0)
msg = f"最小值" + ", ".join([f"{x:.4f}" for x in ret.x])
msg += f"\nf(x)={ret.fun[0]:.4f}"
print(msg)
'''
最小值0.9557, 0.5371, 1.5714, 1.6931, 5.2294
f(x)=0.0000
'''
到了这里，关于【scipy】Python调用非线性最小二乘法的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
