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

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

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

【關鍵詞】SciPy? 最小二乘法? 擬合

1、引言

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

2、理論基礎

2.1最小二乘法理論基礎

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

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

2.2 SciPy基礎

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

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

3、最小二乘法求解設計實現(xiàn)

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

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

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

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

3.1.2實現(xiàn)程序

import numpy as np

from scipy import linalg

import matplotlib.pyplot as plt

import time

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

plt.rcParams['axes.unicode_minus'] = False? ? ? # 用來正常顯示負號

start = time.time（）

np.random.seed（0）

x=np.linspace（-3，3，150）? ? ? ? ? ? ? #x的定義域，150個模擬的未擾動測試樣本

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

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

界干擾