[数值计算-17]：最小二乘法的求解2 - 二元二次线性方程组求解_二次多项式的最小二乘法

作者主页(文火冰糖的硅基工坊)：https://blog.csdn.net/HiWangWenBing

本文网址：https://blog.csdn.net/HiWangWenBing/article/details/119977443

---

目录

第1章 什么是函数逼近？

1.1 抛出问题1：函数插值

1.2  抛出问题2：函数拟合

1.3  什么是回归

第2章 函数拟合

2.1 插值有时候插值并不是最好的选择

2.2 什么函数拟合

 2.2 函数拟合的基本步骤

第3章 最小二乘法的线性方程组解法

3.1 什么是最小二乘法

 3.2 线性拟合函数的最小二乘法的几何意义

3.3 最小二乘算法实现

3.4 最小二乘法求解1：“法向量”方程组法

3.5 最小二乘法求解1：二元线性方程组示例

 第4章 最小二乘法求解2：多元线性方程组解法

4.0 前置条件：

4.1 步骤1：构建样本数据集

4.2 步骤2：构建拟合函数

4.3 步骤3：构建最小二乘的残差函数，并计算最佳参数（最关键步骤）

4.4 步骤4：利用拟合函数进行数据预测

4.5 步骤5：图形展示

---

---

第1章 什么是函数逼近？

1.1 抛出问题1：函数插值

利用有限的样本数据，发现其内在的规律，并用这个规律预测未来新的数据。

（1）单个数据点

• 0次函数通过样本点：唯一确定一个点 y = f(x) = a0
• 1次直线函数通过样本点：可以有无数 y = f(x) = a1x + a0
• 2次抛物线函数通过样本点：可以有无数 y = f(x) = a1x^2 + a1x + a0

（2）2个数据点

• 0次函数通过样本点：无
• 1次直线函数通过样本点：唯一直线 y = f(x) = a1x + a0
• 2次抛物线函数通过样本点：可以有无数 y = f(x) = a1x^2 + a1x + a0

（3）3个数据点

• 0次函数通过样本点：无
• 1次直线函数通过样本点：无
• 2次抛物线函数通过样本点：唯一抛物线y = f(x) = a1x^2 + a1x + a0

问题：

如果有(xn+1,yn+1), (xn,yn)........(x1,y1), (x0,y0)样本点，那么如何选择一个最低次的多项式函数，可以穿过上述样本点？

推测：

对应n+1个点，可以唯一的确定一个一元n次的多项式函数，该多项式函数可以穿越所有n+1个点。

1.2  抛出问题2：函数拟合

如果有n个点，不要求选出的函数穿越所有的点，而是根据这些点构建的轮廓，选择一个更低维度（次数）的函数尽可能的靠近这些样本点呢？

 当函数的次数远远小于样本点的次数是，该如何选择低次的函数?

1.3  什么是回归

回归与拟合是基本相同的概念。线性拟合通常称为线性回归。

第2章 函数拟合

2.1 插值有时候插值并不是最好的选择

（1）简单直线拟合好于多项式拟合

在上图图有，有几十个样本点，如果用多项式插值来进行函数逼近，虽然插值函数在样本点的误差为0，但很显然，需要高次多项式来进行逼近。

如果样本点有几百个，那么多项式的次数高达几百次。

从图中可以大致看出，要逼近上述样本点，只需要一个简单的线性函数即可：y=ax+b，这种函数逼近称为拟合。

（2）抛物线逼近

同理，对上述样本点的逼近，也不需要几十次的多项式进行插值，值需要一个简单的抛物线就可以进行逼近：y =ax^2 + bx + c，，这种函数逼近称为拟合。

（3）三角函数拟合

2.2 什么函数拟合

 

 2.2 函数拟合的基本步骤

步骤1：构建样本

（1）读取已有样本数据，并图形展示样本数据。或者

（2）构建自行测试样本数据

（3）图形化展示样布数据（散点图）

步骤2：选择函数模型

根据散点图，选择函数类，函数类可以从初等函数中进行选取，如线性函数、二次或多次多项式函数、三角函数等。

步骤3： 构建最小二乘的残差函数，并计算最佳参数（最关键步骤）

（1）选择最佳拟合的范数（这里选择最小二乘法）

选取了函数类型后，每个函数都有自己的待定参数，不同的参数，其拟合效果是不同的，如何选择函数的参数，使得拟合效果最好呢？常有的有最小二乘法。

最小二乘法是函数拟合最重要的环节之一。

（2）计算最佳参数（这是最关键、最核心的步骤）

最小二乘法得到的loss残差函数是一个多元二次多项式函数：

• 残差函数的“元”的个数，就是拟合函数的参数
• 残差函数的“元”的次数：2次
• 拟合最好的依据：多元函数的最小值点处的各个变量值，就是拟合函数最好的参数。

求解多元二次多项式函数的方法有多种：

• 法向量解线性方程组法
• 梯度下降法

无论哪种方法，都是通过求残差函数的最小值，基于指定样本数据，获得了最佳参数的拟合函数。

步骤4：利用拟合函数进行数据预测

生成一组输入数据，利用第4步获得的拟合函数的参数以及对应的拟合函数，对数据进行函数运算（预测），得到新的输出数据。

步骤5：图形展示

把预测数据与采样数据、理论曲线一起进行图形展示。

第3章 最小二乘法的线性方程组解法

3.1 什么是最小二乘法

 

 3.2 线性拟合函数的最小二乘法的几何意义

注意： 

最小二乘选择一个适当的直线参数，使得所有样本点处的yi - f(xi) = yi -（axi+b)的平方和的平均值最小，而不是所有样本点到拟合直线的距离的和的平均值最小。

这样选择，主要处于如下原因：

（1）计算简单：求点到拟合曲线的距离要比直接计算样本点处的Y值的差复杂很多。

（2）符合拟合函数原本的物理意义：拟合函数要在所有的样本点xi处，其f(xi）也尽可能的相等，或者说误差尽可能的小。

3.3 最小二乘算法实现

（1）最小二乘法求解1：一元参数的解析法求解

（2）最小二乘法求解2：“法向量”线性方程组求解法

（3）最小二乘法求解3：梯度下降法

最小二乘（二范数）的优点：

 一阶和二阶导数的存在，使得该残差函数，很容易通过“梯度下降法”获得其最小值和对应的各个参数的值。

3.4 最小二乘法求解1：“法向量”方程组法

 

 

 

 

 

3.5 最小二乘法求解1：二元线性方程组示例

 

 第4章 最小二乘法求解2：多元线性方程组解法

4.0 前置条件：

#导入库
from math import *
import time
import numpy as np
import matplotlib.pyplot as plt #画图工具
from pylab import mpl           #中文字体
from scipy import optimize      #最小二乘算法的算法库

4.1 步骤1：构建样本数据集

#步骤1：构建样本
 
#(1) 采用np, 直接手工生成样本的输入：一组等距离的分布在[-1,1]之间的100个点
sample_numbers = 50
 
x_data = np.linspace(0, 1, sample_numbers)
 
#(2) 为这些数据手工打上理论输出值（标签值）：y = 2x + 1
y_data_pure =  2 * x_data + 1.0
 
#（3）为了模拟现实情况，通过随机数来模拟数据噪声
noise_range = 0.4
np.random.seed(10) #设置随机种子, 确保不同时候，执行结果是相同的
#randn(n)生成的0为均值，1为标准差的正态分布的n个随机数。
y_noise = np.random.randn(*x_data.shape) * noise_range  # *x_data.shape：输入样本的维度或个数
 
#（4）人工生成样本的输出：理论值 + 噪声
y_data_noise = y_data_pure + y_noise
 
#(5) 显示样本数据
# 样本的散点图
plt.scatter(x_data, y_data_noise, label="sample", color="black")
 
# 内在的、理论的曲线图
plt.plot(x_data, y_data_pure, label="f_pure(x)", color="blue", linewidth = 4)
 
#设置属性
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.title("线性拟合")
plt.legend(loc="upper left")
plt.show()

4.2 步骤2：构建拟合函数

#步骤2：构建拟合函数：二元一次拟合函数
def f_line_wb(x, w, b):
    return (w*x + b)

4.3 步骤3：构建最小二乘的残差函数，并计算最佳参数（最关键步骤）

#步骤3-1：利用python库提供的最小二乘算法来计算拟合函数的参数
print("使用无噪声数据：")
popt, pcov = optimize.curve_fit(f_line_wb, x_data, y_data_pure)
print(popt)
print(pcov)
w_scipy =  popt[0]
b_scipy =  popt[1]
print("参数w=", w_scipy)
print("参数b=", b_scipy)
 
print("\n使用有噪声数据：")
popt, pcov = optimize.curve_fit(f_line_wb, x_data, y_data_noise)
print(popt)
print(pcov)
w_scipy =  popt[0]
b_scipy =  popt[1]
print("参数w=", w_scipy)
print("参数b=", b_scipy)

使用无噪声数据：
[2. 1.]
[[ 0. -0.]
 [-0.  0.]]
参数w= 2.0
参数b= 1.0

使用有噪声数据：
[1.91826746 1.08186076]
[[ 0.0310567  -0.01552835]
 [-0.01552835  0.01045787]]
参数w= 1.9182674578022025
参数b= 1.0818607577986927

备注说明：

这里利用python库提供的最小二乘算法来计算拟合函数的参数。

因此没有展现用数值的方法求最小二乘最佳参数的过程。

如下函数就展现通过求解多元线性方程组的方法，获得拟合函数的参数a和b。

#步骤3-2： 自定义最小二乘求解拟合函数参数：偏导+线性方程组（求偏导后一次函数）
def usr_curve_fit(input_f, input_x_data, input_y_data_noise):
    #获取采样数据的个数
    n = len(input_x_data)
    
    #计算线性方程组矩阵的数据
    sum_xi = np.sum(input_x_data)
    sum_yi = np.sum(input_y_data_noise)
    sum_xiyi = np.sum(input_x_data * input_y_data_noise)
    sum_xixi = np.sum(input_x_data**2)
    
    #构建线性方程组矩阵
    A=[[0,0],[0,0]]
    b=[0,0]
    A[0][0] = n
    A[0][1] = sum_xi
    A[1][0] = sum_xi
    A[1][1] = sum_xixi
    b[0] = sum_yi
    b[1] = sum_xiyi
    
    #解方程组
    rst = np.linalg.solve(A,b)
    
    #返回拟合函数参数
    return ([rst[1],rst[0]])
 
 
print("使用无噪声数据：")
popt = usr_curve_fit(f_line_wb, x_data, y_data_pure)
print(popt)
print(pcov)
w_usr =  popt[0]
b_usr =  popt[1]
print("参数w=", popt[0])
print("参数b=", popt[1])
 
print("\n使用有噪声数据：")
popt = usr_curve_fit(f_line_wb, x_data, y_data_noise)
w_usr =  popt[0]
b_usr =  popt[1]
print("参数w=", w_usr)
print("参数b=", b_usr)

使用无噪声数据：
[2.0, 0.9999999999999997]
[[ 0.0310567  -0.01552835]
 [-0.01552835  0.01045787]]
参数w= 2.0
参数b= 0.9999999999999997

使用有噪声数据：
参数w= 1.9182674564112352
参数b= 1.081860759861322

可以看出，通过自定义的线性方程组求解获取的拟合参数的参数与通过scipy获得的参数值是一致的。

scipy库：参数a= 1.9182674578022025    参数b= 1.0409944876382622

自定义 ：参数a= 1.9182674564112352    参数b= 1.0818607577986927

4.4 步骤4：利用拟合函数进行数据预测

# 步骤4：利用获得的拟合函数进行数据预测
print("scipy:", w_scipy,b_scipy)
print("user :", w_usr, b_usr)
 
# scipy算法的拟合数据
y_data_scipy = f_line_wb(x_data, w_scipy, b_scipy)
 
# 线性方程组求解的拟合数据
y_data_usr  = f_line_wb(x_data, w_usr, b_usr)

scipy: 1.9182674578022025 1.0818607577986927
user : 1.9182674564112352 1.081860759861322

4.5 步骤5：图形展示

#步骤5： 图形化展示
#(1) 显示样本数据曲线
plt.scatter(x_data, y_data_noise, label="sample", color="black")
 
#(2) 显示理论数据曲线
plt.plot(x_data, y_data_pure, label="intrinsic", color="blue", linewidth = 2)
 
#(3-1) 显示预测数据曲线 -  scipy库实现
plt.plot(x_data, y_data_scipy,  label="predict",  color="red",  linewidth = 2)
 
#(3-2) 显示预测数据曲线 -  自定义实现
plt.plot(x_data, y_data_usr,  label="predict",  color="green",  linewidth = 2)
 
#(3-3) 显示预测数据曲线 -  自定义梯度下降法
#plt.plot(x_data, y_data_grad,  label="predict",  color="green",  linewidth = 2)
 
#设置属性
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
plt.title("线性拟合")
plt.legend(loc="upper left")
plt.show()

 从上图可以看出：线性方程组求解获得拟合函数与scipy库获得的拟合函数，基本重合。

---

作者主页(文火冰糖的硅基工坊)：https://blog.csdn.net/HiWangWenBing

本文网址：https://blog.csdn.net/HiWangWenBing/article/details/119977443