二、 拟合与优化

  1. scipy的optimize模块提供了许多数值优化算法。

  2. 求解非线性方程组：

    1.   scipy.optimize.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0,
    3.   xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)
    • func：是个可调用对象，它代表了非线性方程组。给他传入方程组的各个参数，它返回各个方程残差（残差为零则表示找到了根）
    • x0：预设的方程的根的初始值
    • args：一个元组，用于给func提供额外的参数。
    • fprime：用于计算func的雅可比矩阵（按行排列）。默认情况下，算法自动推算
    • full_output：如果为True，则给出更详细的输出
    • col_deriv：如果为True，则计算雅可比矩阵更快（按列求导）
    • xtol：指定算法收敛的阈值。当误差小于该阈值时，算法停止迭代
    • maxfev：设定算法迭代的最大次数。如果为零，则为 100*(N+1)，N为x0的长度
    • band：If set to a two-sequence containing the number of sub- and super-diagonals within the band of the Jacobi matrix, the Jacobi matrix is considered banded (only for fprime=None)
    • epsfcn：采用前向差分算法求解雅可比矩阵时的步长。
    • factor：它决定了初始的步长
    • diag：它给出了每个变量的缩放因子

    返回值：

    • x：方程组的根组成的数组

    • infodict：给出了可选的输出。它是个字典，其中的键有：

      • nfev：func调用次数
      • njev：雅可比函数调用的次数
      • fvec：最终的func输出
      • fjac：the orthogonal matrix, q, produced by the QR factorization of the final approximate Jacobian matrix, stored column wise
      • r：upper triangular matrix produced by QR factorization of the same matrix

    • ier：一个整数标记。如果为 1，则表示根求解成功

    • mesg：一个字符串。如果根未找到，则它给出详细信息

    假设待求解的方程组为：

    二、 拟合与优化 - 图1

    那么我们的func函数为：

    1.   def func(x):
    3.     x1,x2,x3=x.tolist() # x 为向量，形状为 (3,)
    5.     return np.array([f1(x1,x2,x3),f2(x1,x2,x3),f3(x1,x2,x3)])

    数组的.tolist()方法能获得标准的python列表

    而雅可比矩阵为：

    二、 拟合与优化 - 图2

    1.   def fprime(x):
    3.     x1,x2,x3=x.tolist() # x 为向量，形状为 (3,)
    5.      return np.array([[df1/dx1,df1/dx2,df1/df3],
    7.          [df2/dx1,df2/dx2,df2/df3],
    9.          [df3/dx1,df3/dx2,df3/df3]])

    fsolve

  3. 最小二乘法拟合数据：

    1.   scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0,
    3.   ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, 
    5.   factor=100, diag=None)
    • func：是个可调用对象，给出每次拟合的残差。最开始的参数是待优化参数；后面的参数由args给出
    • x0：初始迭代值
    • args：一个元组，用于给func提供额外的参数。
    • Dfun：用于计算func的雅可比矩阵（按行排列）。默认情况下，算法自动推算。它给出残差的梯度。最开始的参数是待优化参数；后面的参数由args给出
    • full_output：如果非零，则给出更详细的输出
    • col_deriv：如果非零，则计算雅可比矩阵更快（按列求导）
    • ftol：指定相对的均方误差的阈值
    • xtol：指定参数解收敛的阈值
    • gtol：Orthogonality desired between the function vector and the columns of the Jacobian
    • maxfev：设定算法迭代的最大次数。如果为零：如果为提供了Dfun，则为 100*(N+1)，N为x0的长度；如果未提供Dfun,则为200*(N+1)
    • epsfcn：采用前向差分算法求解雅可比矩阵时的步长。
    • factor：它决定了初始的步长
    • diag：它给出了每个变量的缩放因子

    返回值：

    • x：拟合解组成的数组

    • cov_x：Uses the fjac and ipvt optional outputs to construct an estimate of the jacobian around the solution

    • infodict：给出了可选的输出。它是个字典，其中的键有：

      • nfev：func调用次数
      • fvec：最终的func输出
      • fjac：A permutation of the R matrix of a QR factorization of the final approximate Jacobian matrix, stored column wise.
      • ipvt：An integer array of length N which defines a permutation matrix, p, such that fjacp = qr, where r is upper triangular with diagonal elements of nonincreasing magnitude

    • ier：一个整数标记。如果为 1/2/3/4，则表示拟合成功

    • mesg：一个字符串。如果解未找到，则它给出详细信息

    假设我们拟合的函数是 二、 拟合与优化 - 图4，其中 二、 拟合与优化 - 图5 为参数。假设数据点的横坐标为 二、 拟合与优化 - 图6，纵坐标为 二、 拟合与优化 - 图7，那么我们可以给出func为：

    1.   def func(p,x,y):
    3.     a,b,c=p.tolist() # 这里p 为数组，形状为 (3,)； x,y 也是数组，形状都是 (N,)
    5.     return f(x,y;a,b,c))

    其中 args=(X,Y)

    而雅可比矩阵为 二、 拟合与优化 - 图8，给出Dfun为：

    1.   def func(p,x,y):
    3.     a,b,c=p.tolist()
    5.     return np.c_[df/da,df/db,df/dc]# 这里p为数组，形状为 (3,)；x,y 也是数组，形状都是 (N,)

    其中 args=(X,Y)

    leastsq

  4. scipy提供了另一个函数来执行最小二乘法的曲线拟合：

    1.  scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, 
    3.     check_finite=True, bounds=(-inf, inf), method=None, **kwargs)
    • f：可调用函数，它的优化参数被直接传入。其第一个参数一定是xdata，后面的参数是待优化参数
    • xdata：x坐标
    • ydata：y坐标
    • p0：初始迭代值
    • sigma：y值的不确定性的度量
    • absolute_sigma： If False, sigma denotes relative weights of the data points. The returned covariance matrix pcov is based on estimated errors in the data, and is not affected by the overall magnitude of the values in sigma. Only the relative magnitudes of the sigma values matter.If True, sigma describes one standard deviation errors of the input data points. The estimated covariance in pcov is based on these values.
    • check_finite：如果为True，则检测输入中是否有nan或者inf
    • bounds：指定变量的取值范围
    • method：指定求解算法。可以为 'lm'/'trf'/'dogbox'
    • kwargs：传递给 leastsq/least_squares的关键字参数。

    返回值：

    • popt：最优化参数
    • pcov：The estimated covariance of popt.

    假设我们拟合的函数是 二、 拟合与优化 - 图10，其中 二、 拟合与优化 - 图11 为参数。假设数据点的横坐标为 二、 拟合与优化 - 图12，纵坐标为 二、 拟合与优化 - 图13，那么我们可以给出func为：

    1.   def func(x,a,b,c):
    3.     return f(x;a,b,c)#x 为数组，形状为 (N,)

    curve_fit

  5. 求函数最小值：

    1.   scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, 
    3.   hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

    • fun：可调用对象，待优化的函数。最开始的参数是待优化的自变量；后面的参数由args给出

    • x0：自变量的初始迭代值

    • args：一个元组，提供给fun的额外的参数

    • method：一个字符串，指定了最优化算法。可以为：'Nelder-Mead'、'Powell'、'CG'、 'BFGS'、'Newton-CG'、'L-BFGS-B'、'TNC'、'COBYLA'、'SLSQP'、 'dogleg'、'trust-ncg'

    • jac：一个可调用对象（最开始的参数是待优化的自变量；后面的参数由args给出），雅可比矩阵。只在CG/BFGS/Newton-CG/L-BFGS-B/TNC/SLSQP/dogleg/trust-ncg算法中需要。如果jac是个布尔值且为True，则会假设fun会返回梯度；如果是个布尔值且为False，则雅可比矩阵会被自动推断（根据数值插值）。

    • hess/hessp：可调用对象（最开始的参数是待优化的自变量；后面的参数由args给出），海森矩阵。只有Newton-CG/dogleg/trust-ncg算法中需要。二者只需要给出一个就可以，如果给出了hess，则忽略hessp。如果二者都未提供，则海森矩阵自动推断

    • bounds：一个元组的序列，给定了每个自变量的取值范围。如果某个方向不限，则指定为None。每个范围都是一个(min,max)元组。

    • constrants：一个字典或者字典的序列，给出了约束条件。只在COBYLA/SLSQP中使用。字典的键为：

      • type：给出了约束类型。如'eq'代表相等；'ineq'代表不等
      • fun：给出了约束函数
      • jac：给出了约束函数的雅可比矩阵（只用于SLSQP）
      • args：一个序列，给出了传递给fun和jac的额外的参数

    • tol：指定收敛阈值

    • options：一个字典，指定额外的条件。键为：

      • maxiter：一个整数，指定最大迭代次数
      • disp：一个布尔值。如果为True，则打印收敛信息
    • callback：一个可调用对象，用于在每次迭代之后调用。调用参数为x_k，其中x_k为当前的参数向量

    返回值：返回一个OptimizeResult对象。其重要属性为：

    • x：最优解向量
    • success：一个布尔值，表示是否优化成功
    • message：描述了迭代终止的原因

    假设我们要求解最小值的函数为： 二、 拟合与优化 - 图15，则雅可比矩阵为：

    二、 拟合与优化 - 图16

    则海森矩阵为：

    二、 拟合与优化 - 图17

    于是有：

    1.   def fun(p):
    3.     x,y=p.tolist()#p 为数组，形状为 (2,)
    5.     return f(x,y)
    7.   def jac(p):
    9.     x,y=p.tolist()#p 为数组，形状为 (2,)
    11.     return np.array([df/dx,df/dy])
    13.   def hess(p):
    15.     x,y=p.tolist()#p 为数组，形状为 (2,)
    17.     return np.array([[ddf/dxx,ddf/dxdy],[ddf/dydx,ddf/dyy]])

    minimize

  6. 常规的最优化算法很容易陷入局部极值点。basinhopping算法是一个寻找全局最优点的算法。

    1. scipy.optimize.basinhopping(func, x0, niter=100, T=1.0, stepsize=0.5, 
    3.   minimizer_kwargs=None,take_step=None, accept_test=None, callback=None, 
    5.   interval=50, disp=False, niter_success=None)
    • func：可调用函数。为待优化的目标函数。最开始的参数是待优化的自变量；后面的参数由minimizer_kwargs字典给出
    • x0：一个向量，设定迭代的初始值
    • niter：一个整数，指定迭代次数
    • T：一个浮点数，设定了“温度”参数。
    • stepsize：一个浮点数，指定了步长
    • minimizer_kwargs：一个字典，给出了传递给scipy.optimize.minimize的额外的关键字参数。
    • take_step：一个可调用对象，给出了游走策略
    • accept_step：一个可调用对象，用于判断是否接受这一步
    • callback：一个可调用对象，每当有一个极值点找到时，被调用
    • interval：一个整数，指定stepsize被更新的频率
    • disp：一个布尔值，如果为True，则打印状态信息
    • niter_success：一个整数。Stop the run if the global minimum candidate remains the same for this number of iterations.

    返回值：一个OptimizeResult对象。其重要属性为：

    • x：最优解向量
    • success：一个布尔值，表示是否优化成功
    • message：描述了迭代终止的原因

    假设我们要求解最小值的函数为： 二、 拟合与优化 - 图19，于是有：

    1.   def fun(p):
    3.     x,y=p.tolist()#p 为数组，形状为 (2,)
    5.     return f(x,y)

    basinhopping