文章目錄

• • 簡介與構(gòu)造函數(shù)
  • 迭代策略
  • 雅可比矩陣
  • 測試

簡介與構(gòu)造函數(shù)

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

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

其中$\rho$表示損失函數(shù)，可以理解為對$f_i(x)$的一次預(yù)處理。

scipy.optimize中封裝了非線性最小二乘法函數(shù)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為必選參數(shù)，func為待求解函數(shù)，x0為函數(shù)輸入的初值，這兩者無默認(rèn)值，為必須輸入的參數(shù)。

bound為求解區(qū)間，默認(rèn)$(?\infty ,\infty )$，verbose為1時(shí)，會(huì)有終止輸出，為2時(shí)會(huì)print更多的運(yùn)算過程中的信息。此外下面幾個(gè)參數(shù)用于控制誤差，比較簡單。

默認(rèn)值備注

ftol	$10^{?8}$	函數(shù)容忍度
xtol	$10^{?8}$	自變量容忍度
gtol	$10^{?8}$	梯度容忍度
x_scale	1.0	變量的特征尺度
f_scale	1.0	殘差邊際值

loss為損失函數(shù)，就是上面公式中的$\rho$，默認(rèn)為linear，可選值包括

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

迭代策略

上面的公式僅給出了算法的目的，但并未暴露其細(xì)節(jié)。關(guān)于如何找到最小值，則需要確定搜索最小值的方法，method為最小值搜索的方案，共有三種選項(xiàng)，默認(rèn)為trf

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

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

對于優(yōu)化問題${\min }_{x\in R^n}f(x)$而言，定義當(dāng)前點(diǎn)的鄰域

Ω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$為置信半徑，假設(shè)在這個(gè)鄰域內(nèi)，目標(biāo)函數(shù)可以近似為線性或二次函數(shù)，則可通過二次模型得到區(qū)間中的極小值點(diǎn)$s_k$。然后以這個(gè)極小值點(diǎn)為中心，繼續(xù)優(yōu)化信賴域所對應(yīng)的區(qū)間。

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

$$\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$被定義為黑塞矩陣的近似，當(dāng)其表示為雅可比矩陣$J_k^T{J}_k$時(shí)，所得到的迭代方案便是著名的高斯牛頓法，而LM算法在在高斯牛頓的基礎(chǔ)上，添加了一個(gè)阻尼因子，可記作$J_k^T{J}_k+\mu I$。

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

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

雅可比矩陣

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

此外，可以輸入自定義函數(shù)來計(jì)算雅可比矩陣。

測試

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

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 '''

總結(jié)

以上是生活随笔為你收集整理的【scipy】Python调用非线性最小二乘法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。