文章目錄

• • 簡介與構造函數
  • 迭代策略
  • 雅可比矩陣
  • 測試

簡介與構造函數

在scipy中，非線性最小二乘法的目的是找到一組函數，使得誤差函數的平方和最小，可以表示為如下公式

$$\underset{f_i}{\mathrm{argmin}}F(x)=0.5\sum _{i=0}^{m?1}\rho (f_i(x{)}^2),x\in [L,R]$$

其中$\rho$表示損失函數，可以理解為對$f_i(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	$10^{?8}$	函數容忍度
xtol	$10^{?8}$	自變量容忍度
gtol	$10^{?8}$	梯度容忍度
x_scale	1.0	變量的特征尺度
f_scale	1.0	殘差邊際值

loss為損失函數，就是上面公式中的$\rho$，默認為linear，可選值包括

• linear：$\rho (z)=z$，此為標準非線性最小二乘法
• soft_l1： $\rho (z)=2(\sqrt{1+z}?1)$， 相當于L1損失的平滑
• huber： 當$z?1$時，$\rho (z)=z$，否則$\rho (z)=2\sqrt{z}?1$， 表現與soft_l1相似
• cauchy： $\rho (z)=\ln (1+z)$
• arctan：$\rho (z)=\arctan (z)$

迭代策略

上面的公式僅給出了算法的目的，但并未暴露其細節。關于如何找到最小值，則需要確定搜索最小值的方法，method為最小值搜索的方案，共有三種選項，默認為trf

• trf：即Trust Region Reflective，信賴域反射算法
• dogbox：信賴域狗腿算法
• lm：Levenberg-Marquardt算法

這三種方法都是信賴域方法的延申，信賴域的優化思想其實就是從單點的迭代變成了區間的迭代，由于本文的目的是介紹scipy中所封裝好的非線性最小二乘函數，故而僅對其原理做簡略的介紹。

對于優化問題${\min }_{x\in R^n}f(x)$而言，定義當前點的鄰域

Ωk={x∈Rn∣∥x?xk∥?r}\Omega_k=\{x\in R^n\big\vert \Vert x-x_k\Vert\leqslant r\}Ωk?={x∈Rn?∥x?xk?∥?r}

其中$r$為置信半徑，假設在這個鄰域內，目標函數可以近似為線性或二次函數，則可通過二次模型得到區間中的極小值點$s_k$。然后以這個極小值點為中心，繼續優化信賴域所對應的區間。

記$s=x?x_k$, 表示步長；$g_k=?f(x_k),B_k\approx {?}^{2}f(x_k)$分別是函數的一階導數和黑塞矩陣的近似，則對上述問題進行Taylor展開，可以得到一個二階的近似模型

$$\mathrm{argmin}{q}^k(s)=f(x_k)+g_k^{T}s+\frac{1}{2}{s}^T{B}_{k}s,\Vert s\Vert ?r$$

$B_k$被定義為黑塞矩陣的近似，當其表示為雅可比矩陣$J_k^T{J}_k$時，所得到的迭代方案便是著名的高斯牛頓法，而LM算法在在高斯牛頓的基礎上，添加了一個阻尼因子，可記作$J_k^T{J}_k+\mu I$。

狗腿法的鮑威爾提出的一種迭代方案，特點是新定義了下降比，即隨著$s$的不斷前進，目標函數和模型函數都會發生變化，則可定義二者的變化率$r_k=\frac{f(x_k)?f(x_k+s)}{q^k(0)?q^k(s)}$，根據這個值的變化，來調整信賴域半徑$r$的值。

以上就是信賴域方法的基本原理。

雅可比矩陣

在了解了信賴域方法之后，就會明白雅可比矩陣在數值求解時的重要作用，而如何計算雅可比矩陣，則是接下來需要考慮的問題。jac參數為計算雅可比矩陣的方法，主要提供了三種方案，分別是基于兩點的2-point；基于三點的3-point；以及基于復數步長的cs。一般來說，三點的精度高于兩點，但速度也慢一倍。

此外，可以輸入自定義函數來計算雅可比矩陣。

測試

最后，測試一下非線性最小二乘法

import numpy as np from scipy.optimize import least_squaresdef test(xs):_sum = 0.0for i in range(len(xs)):_sum = _sum + (1-np.cos((xs[i]*i)/5)*(i+1))return _sumx0 = 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调用非线性最小二乘法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。