基于SciPy的最小二乘法求解的研究與實(shí)現(xiàn)

劉海靜　劉暢　石慧升　張露予　趙歡

【摘要】本文主要詳細(xì)介紹了如何應(yīng)用Python語言的SciPy庫(kù)中的兩種方法實(shí)現(xiàn)一次函數(shù)模型的最小二乘法擬合優(yōu)化計(jì)算，通過程序分析及結(jié)果比較得出每種方法的效率及精確性。

【關(guān)鍵詞】SciPy? 最小二乘法? 擬合

1、引言

最小二乘法原理最早出現(xiàn)在Legendre在1805年發(fā)表的《計(jì)算彗星軌道的新方法》的附錄中。Legendre認(rèn)為：“賦予誤差的平方和為極小時(shí)，則意味著在這些誤差之間建立了一種均衡性，它阻止了極端情形所施加的過度影響。這非常好地適用于揭示最接近真實(shí)情形的系統(tǒng)狀態(tài)”。

2、理論基礎(chǔ)

2.1最小二乘法理論基礎(chǔ)

最小二乘法主要應(yīng)用于解決連續(xù)或離散函數(shù)模型的最優(yōu)解問題，是測(cè)量工作以及其它科學(xué)工程領(lǐng)域中，應(yīng)用最早也是最廣泛的算法之一。最小二乘法是一種最常見的曲線擬合數(shù)學(xué)優(yōu)化方法，其原理是將復(fù)雜函數(shù)變?yōu)楹?jiǎn)單函數(shù)或者將離散函數(shù)轉(zhuǎn)變?yōu)檫B續(xù)函數(shù)，其意義在于，簡(jiǎn)化函數(shù)或數(shù)學(xué)模型，將原先不易求解和解析困難的問題轉(zhuǎn)化為容易求解和容易析解的問題。其通常的解析過程需要首先在一個(gè)連續(xù)或者離散的原函數(shù)上采樣若干樣本點(diǎn)，接著用若干個(gè)簡(jiǎn)單的函數(shù)與采樣樣本點(diǎn)進(jìn)行擬合，并求擬合函數(shù)與采樣點(diǎn)之間的方差，循環(huán)重復(fù)這一過程，直到樣本中出現(xiàn)一個(gè)最方差，這條擬合曲線就是通過最小二乘法得到的最優(yōu)擬合曲線。

在生產(chǎn)和實(shí)踐過程中，經(jīng)常會(huì)遇到利用一組觀測(cè)數(shù)據(jù)通過最小二乘法來估計(jì)某些未知參數(shù)的問題。最小二乘法在計(jì)算過程中存在一定的誤差，并且使用人工計(jì)算的方法計(jì)算量非常大，但隨著近代計(jì)算機(jī)技術(shù)的不斷發(fā)展，得益于計(jì)算機(jī)強(qiáng)大的計(jì)算能力，最小二乘法計(jì)算量大的問題得到了很好的解決，因此最小二乘法已經(jīng)成為了工程和應(yīng)用數(shù)學(xué)領(lǐng)域最常用的函數(shù)擬合優(yōu)化方法。

2.2 SciPy基礎(chǔ)

SciPy起源于2001年，是一個(gè)面向Python語言的開源科學(xué)計(jì)算數(shù)據(jù)包，常用用于數(shù)學(xué)、科學(xué)、工程領(lǐng)域，可以處理插值、積分、優(yōu)化、整合以及線性代數(shù)模塊、傅里葉變換、圖像處理、常微分方程數(shù)值解的求解、信號(hào)處理等各種數(shù)值計(jì)算的問題。SciPy庫(kù)是構(gòu)建在NumPy庫(kù)的基礎(chǔ)之上的，SciPy提供了許多操作NumPy數(shù)組的函數(shù)，可以有效得計(jì)算NumPy庫(kù)中的矩陣和數(shù)組，使得NumPy庫(kù)和Scipy庫(kù)協(xié)同工作，高效得解決數(shù)據(jù)計(jì)算和數(shù)據(jù)處理問題。

scipy.sparse模塊構(gòu)建了7種稀疏矩陣數(shù)據(jù)結(jié)構(gòu)，7種結(jié)構(gòu)在計(jì)算機(jī)中的以稀疏格式儲(chǔ)存，使得科學(xué)計(jì)算和數(shù)據(jù)處理有效得節(jié)省了空間，提高了運(yùn)行效率。其中最常用的CSR與CSC這兩種方法提供了快速的主軸索引方式和稀疏模式下的快速矩陣向量乘法，大大提高了矩陣的運(yùn)算速度，從而在依托于SciPy庫(kù)的科學(xué)計(jì)算中得到了廣泛的應(yīng)用。scipy.linalg模塊提供了標(biāo)準(zhǔn)的線性代數(shù)運(yùn)算，同時(shí)還加入了NumPy中沒有的高級(jí)函數(shù)，使用BLAS/LAPACK庫(kù)編譯。scipy.optimize模塊提供了非線性方程求解，函數(shù)最小值，曲線擬合和路徑尋找等數(shù)學(xué)優(yōu)化方法。

3、最小二乘法求解設(shè)計(jì)實(shí)現(xiàn)

3.1 Scipy.linalg模擬干擾樣本對(duì)一次函數(shù)求解

本文首先采用scipy.linalg模塊。

3.1.1設(shè)計(jì)思路。利用numpy.linspace在（-3，3）上建立150個(gè)初始數(shù)據(jù)，建立理想模型y=5x+3，通過numpy.random模塊中的隨機(jī)函數(shù)模擬μ為0，σ為0.3的正態(tài)分布外界干擾，并添加種子池，保證每次程序運(yùn)行干擾值固定，對(duì)初始測(cè)試樣本進(jìn)行干擾，通過干擾后的測(cè)試樣本構(gòu)建數(shù)學(xué)模型方程組，抽取方程組中的系數(shù)矩陣，并進(jìn)行numpy縱向堆疊，傳入scipy.linalg模塊lstsq函數(shù)，由于lstsq函數(shù)傳入?yún)?shù)需要為數(shù)組形式，所以對(duì)堆疊之后的結(jié)果進(jìn)行轉(zhuǎn)置，再將轉(zhuǎn)置后的結(jié)果傳入lstsq函數(shù)。

lstsq函數(shù)的返回值共有4個(gè)，所求的結(jié)果系數(shù)在第一個(gè)返回值中，而第一個(gè)返回值為數(shù)組類型，根據(jù)數(shù)學(xué)模型，可以對(duì)應(yīng)相應(yīng)的索引號(hào)打印顯示出最終的最優(yōu)模型系數(shù)，最后通過Matplotlib庫(kù)進(jìn)行數(shù)據(jù)可視化輸出干擾后的測(cè)試樣本和擬合曲線，并在lstsq函數(shù)求解的首尾加入計(jì)時(shí)函數(shù)，具體流程圖如圖1所示。

3.1.2實(shí)現(xiàn)程序

import numpy as np

from scipy import linalg

import matplotlib.pyplot as plt

import time

plt.rcParams['font.sans-serif'] = ['SimHei']? ? ? ?# 用來正常顯示中文標(biāo)簽

plt.rcParams['axes.unicode_minus'] = False? ? ? # 用來正常顯示負(fù)號(hào)

start = time.time（）

np.random.seed（0）

x=np.linspace（-3，3，150）? ? ? ? ? ? ? #x的定義域，150個(gè)模擬的未擾動(dòng)測(cè)試樣本

y=5*x+3? ? ? ? ? ? ? ? ? ? ? ? ? #建立的理想數(shù)學(xué)模型，預(yù)期值

x1=x+np.random.normal（0，0.3，150）? ? #現(xiàn)實(shí)外界條件干擾，通過random函數(shù)模塊模擬外

界干擾