所謂線性最小二乘法，可以理解為是解方程的延續(xù)，區(qū)別在于，當(dāng)未知量遠(yuǎn)小于方程數(shù)的時(shí)候，將得到一個(gè)無(wú)解的問(wèn)題。最小二乘法的實(shí)質(zhì)，是保證誤差最小的情況下對(duì)未知數(shù)進(jìn)行賦值。

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

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

numpy實(shí)現(xiàn)

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

當(dāng)b為M×K的矩陣時(shí)，則對(duì)每一列，都會(huì)計(jì)算一組x。

其返回值共有4個(gè)，分別是擬合得到的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，二者均提供可覆寫(xiě)開(kāi)關(guān)，設(shè)為T(mén)rue可以節(jié)省運(yùn)行時(shí)間，此外，函數(shù)也支持有限性檢查，這是linalg中許多函數(shù)都具備的選項(xiàng)。其返回值與numpy中的最小二乘函數(shù)相同。

cond為浮點(diǎn)型參數(shù)，表示奇異值閾值，當(dāng)奇異值小于cond時(shí)將舍棄。

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

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

速度對(duì)比

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

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

這一次，二者并沒(méi)有拉開(kāi)太大的差距，即使將矩陣維度放大到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

補(bǔ)充

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

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

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

其中ρ表示損失函數(shù)，可以理解為對(duì)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ù)輸入的初值，這兩者無(wú)默認(rèn)值，為必須輸入的參數(shù)。

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

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

迭代策略

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

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

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

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

以上就是信賴(lài)域方法的基本原理。

雅可比矩陣

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

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

測(cè)試

最后，測(cè)試一下非線性最小二乘法

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)用實(shí)現(xiàn)最小二乘法的方法詳解的文章就介紹到這了,更多相關(guān)Python最小二乘法內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論