所謂線性最小平方法,可以理解為是解方程式的延續,差異在於,當未知量遠小於方程式數的時候,將會得到一個無解的問題。最小平方法的實質,是保證誤差最小的情況下對未知數進行賦值。
最小平方法是非常經典的演算法,而且這個名字我們在高中的時候就已經接觸了,屬於極為常用的演算法。先前曾經寫過線性最小平方法的原理,並用Python實現:最小二乘法及其Python實現;以及scipy中非線性最小二乘法的調用方式:非線性最小二乘法(文末補充內容);還有稀疏矩陣的最小平方法:稀疏矩陣最小平方法。
下面講對numpy和scipy中實現的線性最小平方法進行說明,並比較二者的速度。
numpy中便實現了最小二乘法,即lstsq(a,b)用於求解類似於a@x=b中的x,其中,a為M× N的矩陣;則當b為M行的向量時,剛好相當於解線性方程組。對於Ax=b這樣的方程組,如果A是滿秩仿真,那麼可以表示為x=A−1b,否則可以表示為x=(ATA)−1ATb。
當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.linalg同樣提供了最小平方法函數,函數名稱同樣是lstsq,其參數列表為
lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)
其中a, b即Ax= b,二者皆提供可覆寫開關,設為True可以節省運行時間,此外,函數也支援有限性檢查,這是linalg中許多函數都具備的選項。其傳回值與numpy中的最小平方法函數相同。
cond為浮點型參數,表示奇異值閾值,當奇異值小於cond時將捨棄。
lapack_driver為字串選項,表示選用何種LAPACK中的演算法引擎,可選'gelsd', 'gelsy', 'gelss'。
import scipy.linalg as sl xhat1 = sl.lstsq(M, y) print(xhat1[0]) # [0. 1. 2. 3.]
最後,對著兩組最小二乘函數做一個速度上的對比
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呼叫非線性最小平方法
簡介與建構子
在在scipy中,非線性最小二乘法的目的是找到一組函數,使得誤差函數的平方和最小,可以表示為如下公式
其中ρ表示損失函數,可以理解為對fi(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為求解區間,預設(−∞,∞),verbose為1時,會有終止輸出,為2時會print更多的運算過程中的資訊。另外下面幾個參數用來控制誤差,比較簡單。
預設值 | #備註 | |
---|---|---|
ftol | 10-8 | 函數容忍度 |
#xtol | 10-8 | 自變數容忍度 |
gtol | 10-8 | 梯度容忍度 |
x_scale | 1.0 | 變數的特徵尺度 |
1.0 | 殘差邊際值 |