最小二乘法是非常經(jīng)典的算法，而且這個名字我們在高中的時候就已經(jīng)接觸了，屬于極其常用的算法。此前曾經(jīng)寫過線性最小二乘法的原理，并用Python實現(xiàn)：最小二乘法及其Python實現(xiàn)；以及scipy中非線性最小二乘法的調(diào)用方式：非線性最小二乘法（文末補充內(nèi)容）；還有稀疏矩陣的最小二乘法：稀疏矩陣最小二乘法。

下面講對numpy和scipy中實現(xiàn)的線性最小二乘法進行說明，并比較二者的速度。

numpy實現(xiàn)

numpy中便實現(xiàn)了最小二乘法，即lstsq(a,b)用于求解類似于a@x=b中的x，其中，a為M×N的矩陣；則當b為M行的向量時，剛好相當于求解線性方程組。對于Ax=b這樣的方程組，如果A是滿秩仿真，那么可以表示為x=A⁻¹b，否則可以表示為x=(A^TA)⁻¹A^Tb。

當b為M×K的矩陣時，則對每一列，都會計算一組x。

其返回值共有4個，分別是擬合得到的x、擬合誤差、矩陣a的秩、以及矩陣a的單值形式。

import numpy as np
np.random.seed(42)
M = np.random.rand(4,4)
x = np.arange(4)
y = M@x
xhat = np.linalg.lstsq(M,y)
print(xhat[0])
#[0. 1. 2. 3.]

scipy封裝

scipy.linalg同樣提供了最小二乘法函數(shù)，函數(shù)名同樣是lstsq，其參數(shù)列表為

lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)

其中a, b即Ax=b，二者均提供可覆寫開關(guān)，設(shè)為True可以節(jié)省運行時間，此外，函數(shù)也支持有限性檢查，這是linalg中許多函數(shù)都具備的選項。其返回值與numpy中的最小二乘函數(shù)相同。

cond為浮點型參數(shù)，表示奇異值閾值，當奇異值小于cond時將舍棄。

lapack_driver為字符串選項，表示選用何種LAPACK中的算法引擎，可選'gelsd', 'gelsy', 'gelss'。

import scipy.linalg as sl
xhat1 = sl.lstsq(M, y)
print(xhat1[0])
# [0. 1. 2. 3.]

速度對比

最后，對著兩組最小二乘函數(shù)做一個速度上的對比

from timeit import timeit
N = 100
A = np.random.rand(N,N)
b = np.arange(N)

timeit(lambda:np.linalg.lstsq(A, b), number=10)
# 0.015487500000745058
timeit(lambda:sl.lstsq(A, b), number=10)
# 0.011151800004881807

這一次，二者并沒有拉開太大的差距，即使將矩陣維度放大到500，二者也是半斤八兩。

N = 500
A = np.random.rand(N,N)
b = np.arange(N)

timeit(lambda:np.linalg.lstsq(A, b), number=10)
0.389679799991427
timeit(lambda:sl.lstsq(A, b), number=10)
0.35642060000100173

補充

Python調(diào)用非線性最小二乘法

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

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

其中ρ表示損失函數(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ù)輸入的初值，這兩者無默認值，為必須輸入的參數(shù)。

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

	默認值	備注
ftol	10^-8	函數(shù)容忍度
xtol	10^-8	自變量容忍度
gtol	10^-8	梯度容忍度
x_scale	1.0	變量的特征尺度
f_scale	1.0	殘差邊際值

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

迭代策略

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

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

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

其中r為置信半徑，假設(shè)在這個鄰域內(nèi)，目標函數(shù)可以近似為線性或二次函數(shù)，則可通過二次模型得到區(qū)間中的極小值點s_k。然后以這個極小值點為中心，繼續(xù)優(yōu)化信賴域所對應(yīng)的區(qū)間。

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

雅可比矩陣

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

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

測試

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

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

到此這篇關(guān)于Python調(diào)用實現(xiàn)最小二乘法的方法詳解的文章就介紹到這了,更多相關(guān)Python最小二乘法內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: