基于SciPy的最小二乘法求解的研究与实现

刘海静　刘畅　石慧升　张露予　赵欢

【摘要】本文主要详细介绍了如何应用Python语言的SciPy库中的两种方法实现一次函数模型的最小二乘法拟合优化计算，通过程序分析及结果比较得出每种方法的效率及精确性。

【关键词】SciPy  最小二乘法  拟合

1、引言

最小二乘法原理最早出现在Legendre在1805年发表的《计算彗星轨道的新方法》的附录中。Legendre认为：“赋予误差的平方和为极小时，则意味着在这些误差之间建立了一种均衡性，它阻止了极端情形所施加的过度影响。这非常好地适用于揭示最接近真实情形的系统状态”。

2、理论基础

2.1最小二乘法理论基础

最小二乘法主要应用于解决连续或离散函数模型的最优解问题，是测量工作以及其它科学工程领域中，应用最早也是最广泛的算法之一。最小二乘法是一种最常见的曲线拟合数学优化方法，其原理是将复杂函数变为简单函数或者将离散函数转变为连续函数，其意义在于，简化函数或数学模型，将原先不易求解和解析困难的问题转化为容易求解和容易析解的问题。其通常的解析过程需要首先在一个连续或者离散的原函数上采样若干样本点，接着用若干个简单的函数与采样样本点进行拟合，并求拟合函数与采样点之间的方差，循环重复这一过程，直到样本中出现一个最方差，这条拟合曲线就是通过最小二乘法得到的最优拟合曲线。

在生产和实践过程中，经常会遇到利用一组观测数据通过最小二乘法来估计某些未知参数的问题。最小二乘法在计算过程中存在一定的误差，并且使用人工计算的方法计算量非常大，但随着近代计算机技术的不断发展，得益于计算机强大的计算能力，最小二乘法计算量大的问题得到了很好的解决，因此最小二乘法已经成为了工程和应用数学领域最常用的函数拟合优化方法。

2.2 SciPy基础

SciPy起源于2001年，是一个面向Python语言的开源科学计算数据包，常用用于数学、科学、工程领域，可以处理插值、积分、优化、整合以及线性代数模块、傅里叶变换、图像处理、常微分方程数值解的求解、信号处理等各种数值计算的问题。SciPy库是构建在NumPy库的基础之上的，SciPy提供了许多操作NumPy数组的函数，可以有效得计算NumPy库中的矩阵和数组，使得NumPy库和Scipy库协同工作，高效得解决数据计算和数据处理问题。

scipy.sparse模块构建了7种稀疏矩阵数据结构，7种结构在计算机中的以稀疏格式储存，使得科学计算和数据处理有效得节省了空间，提高了运行效率。其中最常用的CSR与CSC这两种方法提供了快速的主轴索引方式和稀疏模式下的快速矩阵向量乘法，大大提高了矩阵的运算速度，从而在依托于SciPy库的科学计算中得到了广泛的应用。scipy.linalg模块提供了标准的线性代数运算，同时还加入了NumPy中没有的高级函数，使用BLAS/LAPACK库编译。scipy.optimize模块提供了非线性方程求解，函数最小值，曲线拟合和路径寻找等数学优化方法。

3、最小二乘法求解设计实现

3.1 Scipy.linalg模拟干扰样本对一次函数求解

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

3.1.1设计思路。利用numpy.linspace在（-3，3）上建立150个初始数据，建立理想模型y=5x+3，通过numpy.random模块中的随机函数模拟μ为0，σ为0.3的正态分布外界干扰，并添加种子池，保证每次程序运行干扰值固定，对初始测试样本进行干扰，通过干扰后的测试样本构建数学模型方程组，抽取方程组中的系数矩阵，并进行numpy纵向堆叠，传入scipy.linalg模块lstsq函数，由于lstsq函数传入参数需要为数组形式，所以对堆叠之后的结果进行转置，再将转置后的结果传入lstsq函数。

lstsq函数的返回值共有4个，所求的结果系数在第一个返回值中，而第一个返回值为数组类型，根据数学模型，可以对应相应的索引号打印显示出最终的最优模型系数，最后通过Matplotlib库进行数据可视化输出干扰后的测试样本和拟合曲线，并在lstsq函数求解的首尾加入计时函数，具体流程图如图1所示。

3.1.2实现程序

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                          #建立的理想数学模型，预期值

x1=x+np.random.normal（0，0.3，150）    #现实外界条件干扰，通过random函数模块模拟外

界干扰